[yt-svn] commit/yt: 14 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Sep 21 20:32:50 PDT 2016
14 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/190c33cab578/
Changeset: 190c33cab578
Branch: yt
User: al007
Date: 2016-09-09 16:31:14+00:00
Summary: Add second order triangle shape functions to mesh_types.yaml
Affected #: 1 file
diff -r c07d12bf5903b40c0c6e65a8ea8f28a82f142e1d -r 190c33cab578d4dcbab5b65a55594b1ee54c326e yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -12,7 +12,7 @@
(1 + x[0])*(1 - x[1])*(1 + x[2])/8.,
(1 + x[0])*(1 + x[1])*(1 + x[2])/8.,
(1 - x[0])*(1 + x[1])*(1 + x[2])/8.]
-
+
Quad4:
mesh_type: Q1
num_dim: 2
@@ -35,4 +35,17 @@
x[1] * (1 - x[2]) / 2.,
(1 - x[0] - x[1]) * (1 + x[2]) / 2.,
x[0] * (1 + x[2]) / 2.,
- x[1] * (1 + x[2]) / 2.]
\ No newline at end of file
+ x[1] * (1 + x[2]) / 2.]
+
+Tri6:
+ mesh_type: T2
+ num_dim: 2
+ num_vertices: 6
+ num_mapped_coords: 2
+ shape_functions: |
+ [1 - 3 * x[0] + 2 * x[0]**2 - 3 * x[1] + 2 * x[1]**2 + 4 * x[0] * x[1],
+ -x[0] + 2 * x[0]**2,
+ -x[1] + 2 * x[1]**2,
+ 4 * x[0] - 4 * x[0]**2 - 4 * x[0] * x[1],
+ 4 * x[0] * x[1],
+ 4 * x[1] - 4 * x[1]**2 - 4 * x[0] * x[1]]
https://bitbucket.org/yt_analysis/yt/commits/a3ad16bcb48e/
Changeset: a3ad16bcb48e
Branch: yt
User: al007
Date: 2016-09-09 17:20:20+00:00
Summary: Add T2Sampler2D to element mappings; also autogenerated element samplers.
Affected #: 3 files
diff -r 190c33cab578d4dcbab5b65a55594b1ee54c326e -r a3ad16bcb48ebb152b89fcba80fad2f219c7ba06 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 T2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void T2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
cdef void W1Function3D(double* fx,
double* x,
double* vertices,
diff -r 190c33cab578d4dcbab5b65a55594b1ee54c326e -r a3ad16bcb48ebb152b89fcba80fad2f219c7ba06 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 T2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = (-x[0] + 2*pow(x[0], 2))*vertices[2] + (-x[1] + 2*pow(x[1], 2))*vertices[4] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[10] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[6] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[0] - phys_x[0] + 4*x[0]*x[1]*vertices[8];
+ fx[1] = (-x[0] + 2*pow(x[0], 2))*vertices[3] + (-x[1] + 2*pow(x[1], 2))*vertices[5] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[11] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[7] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[1] - phys_x[1] + 4*x[0]*x[1]*vertices[9];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void T2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = (-1 + 4*x[0])*vertices[2] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 8*x[0] - 4*x[1])*vertices[6] + 4*x[1]*vertices[8] - 4*x[1]*vertices[10];
+ scol[0] = (-1 + 4*x[1])*vertices[4] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 4*x[0] - 8*x[1])*vertices[10] - 4*x[0]*vertices[6] + 4*x[0]*vertices[8];
+ rcol[1] = (-1 + 4*x[0])*vertices[3] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 8*x[0] - 4*x[1])*vertices[7] + 4*x[1]*vertices[9] - 4*x[1]*vertices[11];
+ scol[1] = (-1 + 4*x[1])*vertices[5] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 4*x[0] - 8*x[1])*vertices[11] - 4*x[0]*vertices[7] + 4*x[0]*vertices[9];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
cdef void W1Function3D(double* fx,
double* x,
double* vertices,
diff -r 190c33cab578d4dcbab5b65a55594b1ee54c326e -r a3ad16bcb48ebb152b89fcba80fad2f219c7ba06 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -1,6 +1,6 @@
"""
This file contains coordinate mappings between physical coordinates and those
-defined on unit elements, as well as doing the corresponding intracell
+defined on unit elements, as well as doing the corresponding intracell
interpolation on finite element data.
@@ -33,8 +33,8 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef double determinant_3x3(double* col0,
- double* col1,
+cdef double determinant_3x3(double* col0,
+ double* col1,
double* col2) nogil:
return col0[0]*col1[1]*col2[2] - col0[0]*col1[2]*col2[1] - \
col0[1]*col1[0]*col2[2] + col0[1]*col1[2]*col2[0] + \
@@ -49,7 +49,7 @@
cdef int i
err = fabs(f[0])
for i in range(1, dim):
- err = fmax(err, fabs(f[i]))
+ err = fmax(err, fabs(f[i]))
return err
@@ -58,7 +58,7 @@
This is a base class for sampling the value of a finite element solution
at an arbitrary point inside a mesh element. In general, this will be done
- by transforming the requested physical coordinate into a mapped coordinate
+ by transforming the requested physical coordinate into a mapped coordinate
system, sampling the solution in mapped coordinates, and returning the result.
This is not to be used directly; use one of the subclasses instead.
@@ -71,11 +71,11 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil:
pass
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -127,43 +127,43 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void map_real_to_unit(self, double* mapped_x,
+ cdef void map_real_to_unit(self, double* mapped_x,
double* vertices, double* physical_x) nogil:
-
+
cdef double[3] col0
cdef double[3] col1
cdef double[3] col2
-
+
col0[0] = vertices[0]
col0[1] = vertices[1]
col0[2] = 1.0
-
+
col1[0] = vertices[2]
col1[1] = vertices[3]
col1[2] = 1.0
-
+
col2[0] = vertices[4]
col2[1] = vertices[5]
col2[2] = 1.0
-
+
det = determinant_3x3(col0, col1, col2)
-
+
mapped_x[0] = ((vertices[3] - vertices[5])*physical_x[0] + \
(vertices[4] - vertices[2])*physical_x[1] + \
(vertices[2]*vertices[5] - vertices[4]*vertices[3])) / det
-
+
mapped_x[1] = ((vertices[5] - vertices[1])*physical_x[0] + \
(vertices[0] - vertices[4])*physical_x[1] + \
(vertices[4]*vertices[1] - vertices[0]*vertices[5])) / det
-
+
mapped_x[2] = 1.0 - mapped_x[1] - mapped_x[0]
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef double sample_at_unit_point(self,
- double* coord,
+ double* coord,
double* vals) nogil:
return vals[0]*coord[0] + vals[1]*coord[1] + vals[2]*coord[2]
@@ -197,16 +197,16 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void map_real_to_unit(self, double* mapped_x,
+ cdef void map_real_to_unit(self, double* mapped_x,
double* vertices, double* physical_x) nogil:
-
+
cdef int i
cdef double d
cdef double[3] bvec
cdef double[3] col0
cdef double[3] col1
cdef double[3] col2
-
+
# here, we express positions relative to the 4th element,
# which is selected by vertices[9]
for i in range(3):
@@ -214,7 +214,7 @@
col0[i] = vertices[0 + i] - vertices[9 + i]
col1[i] = vertices[3 + i] - vertices[9 + i]
col2[i] = vertices[6 + i] - vertices[9 + i]
-
+
d = determinant_3x3(col0, col1, col2)
mapped_x[0] = determinant_3x3(bvec, col1, col2)/d
mapped_x[1] = determinant_3x3(col0, bvec, col2)/d
@@ -225,7 +225,7 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef double sample_at_unit_point(self,
- double* coord,
+ double* coord,
double* vals) nogil:
return vals[0]*coord[0] + vals[1]*coord[1] + \
vals[2]*coord[2] + vals[3]*coord[3]
@@ -265,11 +265,11 @@
u = mapped_coord[1]
v = mapped_coord[2]
w = mapped_coord[0]
- if ((u < thresh) or
- (v < thresh) or
+ if ((u < thresh) or
+ (v < thresh) or
(w < thresh) or
- (fabs(u - 1) < thresh) or
- (fabs(v - 1) < thresh) or
+ (fabs(u - 1) < thresh) or
+ (fabs(v - 1) < thresh) or
(fabs(w - 1) < thresh)):
return 1
return -1
@@ -281,7 +281,7 @@
This is a base class for handling element samplers that require
a nonlinear solve to invert the mapping between coordinate systems.
- To do this, we perform Newton-Raphson iteration using a specified
+ To do this, we perform Newton-Raphson iteration using a specified
system of equations with an analytic Jacobian matrix. This solver
is hard-coded for 3D for reasons of efficiency. This is not to be
used directly, use one of the subclasses instead.
@@ -302,7 +302,7 @@
double* physical_x) nogil:
cdef int i
cdef double d
- cdef double[3] f
+ cdef double[3] f
cdef double[3] r
cdef double[3] s
cdef double[3] t
@@ -313,11 +313,11 @@
# initial guess
for i in range(3):
x[i] = 0.0
-
+
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 3)
-
+ err = maxnorm(f, 3)
+
# begin Newton iteration
while (err > self.tolerance and iterations < self.max_iter):
self.jac(r, s, t, x, vertices, physical_x)
@@ -325,7 +325,7 @@
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)
+ self.func(f, x, vertices, physical_x)
err = maxnorm(f, 3)
iterations += 1
@@ -340,7 +340,7 @@
cdef class Q1Sampler3D(NonlinearSolveSampler3D):
- '''
+ '''
This implements sampling inside a 3D, linear, hexahedral mesh element.
@@ -358,14 +358,14 @@
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
cdef double F, rm, rp, sm, sp, tm, tp
-
+
rm = 1.0 - coord[0]
rp = 1.0 + coord[0]
sm = 1.0 - coord[1]
sp = 1.0 + coord[1]
tm = 1.0 - coord[2]
tp = 1.0 + coord[2]
-
+
F = vals[0]*rm*sm*tm + vals[1]*rp*sm*tm + vals[2]*rp*sp*tm + vals[3]*rm*sp*tm + \
vals[4]*rm*sm*tp + vals[5]*rp*sm*tp + vals[6]*rp*sp*tp + vals[7]*rm*sp*tp
return 0.125*F
@@ -374,10 +374,10 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int check_inside(self, double* mapped_coord) nogil:
- # for hexes, the mapped coordinates all go from -1 to 1
+ # for hexes, the mapped coordinates all go from -1 to 1
# if we are inside the element.
if (fabs(mapped_coord[0]) - 1.0 > self.inclusion_tol or
- fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
+ fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
fabs(mapped_coord[2]) - 1.0 > self.inclusion_tol):
return 0
return 1
@@ -401,7 +401,7 @@
cdef class S2Sampler3D(NonlinearSolveSampler3D):
- '''
+ '''
This implements sampling inside a 3D, 20-node hexahedral mesh element.
@@ -459,7 +459,7 @@
@cython.cdivision(True)
cdef int check_inside(self, double* mapped_coord) nogil:
if (fabs(mapped_coord[0]) - 1.0 > self.inclusion_tol or
- fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
+ fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
fabs(mapped_coord[2]) - 1.0 > self.inclusion_tol):
return 0
return 1
@@ -485,8 +485,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef inline void S2Function3D(double* fx,
- double* x,
- double* vertices,
+ double* x,
+ double* vertices,
double* phys_x) nogil:
cdef int i
cdef double r, s, t, rm, rp, sm, sp, tm, tp
@@ -533,9 +533,9 @@
cdef inline void S2Jacobian3D(double* rcol,
double* scol,
double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
cdef int i
cdef double r, s, t, rm, rp, sm, sp, tm, tp
@@ -616,7 +616,7 @@
cdef class W1Sampler3D(NonlinearSolveSampler3D):
- '''
+ '''
This implements sampling inside a 3D, linear, wedge mesh element.
@@ -648,15 +648,15 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int check_inside(self, double* mapped_coord) nogil:
- # for wedges the bounds of the mapped coordinates are:
+ # for wedges the bounds of the mapped coordinates are:
# 0 <= mapped_coord[0] <= 1 - mapped_coord[1]
# 0 <= mapped_coord[1]
# -1 <= mapped_coord[2] <= 1
if (mapped_coord[0] < -self.inclusion_tol or
mapped_coord[0] + mapped_coord[1] - 1.0 > self.inclusion_tol):
- return 0
+ return 0
if (mapped_coord[1] < -self.inclusion_tol):
- return 0
+ return 0
if (fabs(mapped_coord[2]) - 1.0 > self.inclusion_tol):
return 0
return 1
@@ -675,7 +675,7 @@
near_edge_r = (r < thresh) or (fabs(r + s - 1.0) < thresh)
near_edge_s = (s < thresh)
near_edge_t = fabs(fabs(mapped_coord[2]) - 1.0) < thresh
-
+
# we use ray.instID to pass back whether the ray is near the
# element boundary or not (used to annotate mesh lines)
if (near_edge_r and near_edge_s):
@@ -694,7 +694,7 @@
This is a base class for handling element samplers that require
a nonlinear solve to invert the mapping between coordinate systems.
- To do this, we perform Newton-Raphson iteration using a specified
+ To do this, we perform Newton-Raphson iteration using a specified
system of equations with an analytic Jacobian matrix. This solver
is hard-coded for 2D for reasons of efficiency. This is not to be
used directly, use one of the subclasses instead.
@@ -724,20 +724,20 @@
# initial guess
for i in range(2):
x[i] = 0.0
-
+
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 2)
-
+ err = maxnorm(f, 2)
+
# begin Newton iteration
while (err > 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
- self.func(f, x, vertices, physical_x)
+ self.func(f, x, vertices, physical_x)
err = maxnorm(f, 2)
iterations += 1
@@ -752,7 +752,7 @@
cdef class Q1Sampler2D(NonlinearSolveSampler2D):
- '''
+ '''
This implements sampling inside a 2D, linear, quadrilateral mesh element.
@@ -770,12 +770,12 @@
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
cdef double F, rm, rp, sm, sp
-
+
rm = 1.0 - coord[0]
rp = 1.0 + coord[0]
sm = 1.0 - coord[1]
sp = 1.0 + coord[1]
-
+
F = vals[0]*rm*sm + vals[1]*rp*sm + vals[2]*rp*sp + vals[3]*rm*sp
return 0.25*F
@@ -790,6 +790,49 @@
return 0
return 1
+cdef class T2Sampler2D(NonlinearSolveSampler2D):
+
+ '''
+
+ This implements sampling inside a 2D, linear, quadrilateral mesh element.
+
+ '''
+
+ def __init__(self):
+ super(T2Sampler2D, self).__init__()
+ self.num_mapped_coords = 2
+ self.dim = 2
+ self.func = T2Function2D
+ self.jac = T2Jacobian2D
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
+ cdef double phi0, phi1, phi2, phi3, phi4, phi5
+
+ phi0 = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
+ 2 * coord[1]**2 + 4 * coord[0] * coord[1]
+ phi1 = -coord[0] + 2 * coord[0]**2
+ phi2 = -coord[1] + 2 * coord[1]**2
+ phi3 = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1]
+ phi4 = 4 * coord[0] * coord[1]
+ phi5 = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1]
+
+ return vals[0]*phi0 + vals[1]*phi1 + vals[2]*phi2 + vals[3]*phi3 + \
+ vals[4]*phi4 + vals[5]*phi5
+
+ @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
+
@cython.boundscheck(False)
@cython.wraparound(False)
https://bitbucket.org/yt_analysis/yt/commits/e033d4df50c4/
Changeset: e033d4df50c4
Branch: yt
User: al007
Date: 2016-09-09 22:29:42+00:00
Summary: Edit `check_inside` function for T2Sampler2D; also add some doco
Affected #: 1 file
diff -r a3ad16bcb48ebb152b89fcba80fad2f219c7ba06 -r e033d4df50c462067cecfe593e1ad43bed1cc0fc yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -115,7 +115,8 @@
'''
This implements sampling inside a linear, triangular mesh element.
- This mapping is linear and can be inverted easily.
+ This mapping is linear and can be inverted easily. Note that this
+ implementation uses triangular (or barycentric) coordinates.
'''
@@ -794,7 +795,8 @@
'''
- This implements sampling inside a 2D, linear, quadrilateral mesh element.
+ This implements sampling inside a 2D, quadratic, triangular mesh
+ element. Note that this implementation uses canonical coordinates.
'''
@@ -828,12 +830,13 @@
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
+ cdef int i
+ for i in range(2):
+ if (mapped_coord[i] < -self.inclusion_tol or
+ mapped_coord[i] - 1.0 > self.inclusion_tol):
+ return 0
return 1
-
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
https://bitbucket.org/yt_analysis/yt/commits/7ad10b908654/
Changeset: 7ad10b908654
Branch: yt
User: al007
Date: 2016-09-10 00:49:27+00:00
Summary: Finish implementing second order triange shape function support. Added test.
Affected #: 4 files
diff -r e033d4df50c462067cecfe593e1ad43bed1cc0fc -r 7ad10b90865439f3d029a78b6bb86d92712944dd yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -5,7 +5,7 @@
cimport cython
-
+from libc.math cimport pow
@cython.boundscheck(False)
diff -r e033d4df50c462067cecfe593e1ad43bed1cc0fc -r 7ad10b90865439f3d029a78b6bb86d92712944dd yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -25,7 +25,9 @@
Q1Function2D, \
Q1Jacobian2D, \
W1Function3D, \
- W1Jacobian3D
+ W1Jacobian3D, \
+ T2Function2D, \
+ T2Jacobian2D
cdef extern from "platform_dep.h":
double fmax(double x, double y) nogil
@@ -940,3 +942,20 @@
<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):
+
+ cdef double val
+
+ sampler = T2Sampler2D()
+
+ val = sampler.sample_at_real_point(<double*> vertices.data,
+ <double*> field_values.data,
+ <double*> physical_x.data)
+
+ return val
diff -r e033d4df50c462067cecfe593e1ad43bed1cc0fc -r 7ad10b90865439f3d029a78b6bb86d92712944dd 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
@@ -22,7 +22,8 @@
test_tri_sampler, \
test_quad_sampler, \
test_hex20_sampler, \
- test_wedge_sampler
+ test_wedge_sampler, \
+ test_tri2_sampler
def check_all_vertices(sampler, vertices, field_values):
@@ -127,3 +128,16 @@
field_values = np.array([1., 2., 3., 4., 5., 6.])
check_all_vertices(test_wedge_sampler, vertices, field_values)
+
+def test_T2Sampler2D():
+
+ vertices = np.array([[0.1 , 0.2 ],
+ [0.3 , 0.5 ],
+ [0.2 , 0.9 ],
+ [0.2 , 0.35],
+ [0.25, 0.7 ],
+ [0.15, 0.55]])
+
+ field_values = np.array([15., 37., 49., 32., 46., 24.])
+
+ check_all_vertices(test_tri2_sampler, vertices, field_values)
diff -r e033d4df50c462067cecfe593e1ad43bed1cc0fc -r 7ad10b90865439f3d029a78b6bb86d92712944dd yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -5,9 +5,9 @@
Usage (from the yt/utilities directory):
-python mesh_code_generation.py
+python mesh_code_generation.py
-This will generate the necessary functions and write them to
+This will generate the necessary functions and write them to
yt/utilities/lib/autogenerated_element_samplers.pyx.
@@ -120,7 +120,7 @@
self.jacobian_name = '%sJacobian%dD' % (self.mesh_type, self.num_dim)
if (self.num_dim == 3):
- self.jacobian_header = jac_def_template_3D % self.jacobian_name
+ self.jacobian_header = jac_def_template_3D % self.jacobian_name
self.jacobian_declaration = jac_dec_template_3D % self.jacobian_name
elif (self.num_dim == 2):
@@ -137,15 +137,15 @@
function_code = self.function_header
for i in range(self.num_dim):
function_code += '\t' + ccode(self.f[i, 0], self.fx[i, 0]) + '\n'
-
+
jacobian_code = self.jacobian_header
for i in range(self.num_dim):
jacobian_code += '\t' + ccode(self.J[i,0], self.rcol[i, 0]) + '\n'
jacobian_code += '\t' + ccode(self.J[i,1], self.scol[i, 0]) + '\n'
- if self.num_dim == 2:
+ if self.num_dim == 2:
continue
jacobian_code += '\t' + ccode(self.J[i,2], self.tcol[i, 0]) + '\n'
-
+
return function_code, jacobian_code
def get_interpolator_declaration(self):
@@ -169,9 +169,10 @@
pyx_file.write(file_header)
pyx_file.write("\n \n")
- pyx_file.write("cimport cython \n \n")
+ pyx_file.write("cimport cython \n")
+ pyx_file.write("from libc.math cimport pow \n")
pyx_file.write("\n \n")
-
+
for _, mesh_data in sorted(mesh_types.items()):
codegen = MeshCodeGenerator(mesh_data)
https://bitbucket.org/yt_analysis/yt/commits/f9843e2a1748/
Changeset: f9843e2a1748
Branch: yt
User: al007
Date: 2016-09-12 19:09:41+00:00
Summary: Second order triangle sampler routines for exodus viewing.
- Required additions to element_mappings.pxd and pixelization_routines.pyx.
- Also added @atmyer's line changes fixing vertex fields for non-Hex8 mesh
types in geometry_handler.py.
Affected #: 4 files
diff -r 7ad10b90865439f3d029a78b6bb86d92712944dd -r f9843e2a1748165775d26c9e458f88ff146a4424 yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -172,4 +172,3 @@
@property
def period(self):
return self.ds.domain_width
-
diff -r 7ad10b90865439f3d029a78b6bb86d92712944dd -r f9843e2a1748165775d26c9e458f88ff146a4424 yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -270,7 +270,7 @@
else:
tr = func(self)
if self._cache:
-
+
setattr(self, n, tr)
return tr
return property(cached_func)
@@ -397,7 +397,9 @@
@cached_property
def fcoords_vertex(self):
- ci = np.empty((self.data_size, 8, 3), dtype='float64')
+ nodes_per_elem = self.dobj.index.meshes[0].connectivity_indices.shape[1]
+ dim = self.dobj.ds.dimensionality
+ ci = np.empty((self.data_size, nodes_per_elem, dim), dtype='float64')
ci = YTArray(ci, input_units = "code_length",
registry = self.dobj.ds.unit_registry)
if self.data_size == 0: return ci
@@ -426,10 +428,10 @@
def __iter__(self):
return self
-
+
def __next__(self):
return self.next()
-
+
def next(self):
if len(self.queue) == 0:
for i in range(self.max_length):
@@ -445,4 +447,3 @@
g = self.queue.pop(0)
g._initialize_cache(self.cache.pop(g.id, {}))
return g
-
diff -r 7ad10b90865439f3d029a78b6bb86d92712944dd -r f9843e2a1748165775d26c9e458f88ff146a4424 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -13,7 +13,7 @@
cdef int num_mapped_coords
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -21,7 +21,7 @@
cdef double sample_at_unit_point(self,
double* coord,
double* vals) nogil
-
+
cdef double sample_at_real_point(self,
double* vertices,
@@ -36,7 +36,7 @@
cdef class P1Sampler2D(ElementSampler):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -51,7 +51,7 @@
cdef class P1Sampler3D(ElementSampler):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -67,7 +67,7 @@
# This typedef defines a function pointer that defines the system
# of equations that will be solved by the NonlinearSolveSamplers.
-#
+#
# inputs:
# x - pointer to the mapped coordinate
# vertices - pointer to the element vertices
@@ -78,15 +78,15 @@
# fx - the result of solving the system, should be close to 0
# once it is converged.
#
-ctypedef void (*func_type)(double* fx,
- double* x,
- double* vertices,
+ctypedef void (*func_type)(double* fx,
+ double* x,
+ double* vertices,
double* phys_x) nogil
# This typedef defines a function pointer that defines the Jacobian
-# matrix used by the NonlinearSolveSampler3D. Subclasses needed to
+# matrix used by the NonlinearSolveSampler3D. Subclasses needed to
# define a Jacobian function in this form.
-#
+#
# inputs:
# x - pointer to the mapped coordinate
# vertices - pointer to the element vertices
@@ -98,18 +98,18 @@
# scol - the second column of the jacobian
# tcol - the third column of the jaocobian
#
-ctypedef void (*jac_type3D)(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
+ctypedef void (*jac_type3D)(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
double* phys_x) nogil
# This typedef defines a function pointer that defines the Jacobian
-# matrix used by the NonlinearSolveSampler2D. Subclasses needed to
+# matrix used by the NonlinearSolveSampler2D. Subclasses needed to
# define a Jacobian function in this form.
-#
+#
# inputs:
# x - pointer to the mapped coordinate
# vertices - pointer to the element vertices
@@ -132,19 +132,19 @@
cdef int dim
cdef int max_iter
cdef np.float64_t tolerance
- cdef func_type func
+ cdef func_type func
cdef jac_type3D jac
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
-
+
cdef class Q1Sampler3D(NonlinearSolveSampler3D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -161,7 +161,7 @@
cdef class W1Sampler3D(NonlinearSolveSampler3D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -179,7 +179,7 @@
cdef class S2Sampler3D(NonlinearSolveSampler3D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -199,19 +199,19 @@
cdef int dim
cdef int max_iter
cdef np.float64_t tolerance
- cdef func_type func
+ cdef func_type func
cdef jac_type2D jac
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
-
+
cdef class Q1Sampler2D(NonlinearSolveSampler2D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -221,3 +221,17 @@
double* vals) nogil
cdef int check_inside(self, double* mapped_coord) nogil
+
+cdef class T2Sampler2D(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
diff -r 7ad10b90865439f3d029a78b6bb86d92712944dd -r f9843e2a1748165775d26c9e458f88ff146a4424 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -30,7 +30,8 @@
S2Sampler3D, \
P1Sampler2D, \
Q1Sampler2D, \
- W1Sampler3D
+ W1Sampler3D, \
+ T2Sampler2D
cdef extern from "pixelization_constants.h":
enum:
@@ -201,7 +202,7 @@
my_array[j,i] += (dsp * overlap1) * overlap2
else:
my_array[j,i] = dsp
-
+
return my_array
@cython.cdivision(True)
@@ -309,10 +310,10 @@
cdef np.float64_t r_i, theta_i, dr_i, dtheta_i, dthetamin
cdef np.float64_t costheta, sintheta
cdef int i, pi, pj
-
+
imax = radius.argmax()
rmax = radius[imax] + dradius[imax]
-
+
if input_img is None:
img = np.zeros((buff_size[0], buff_size[1]))
img[:] = np.nan
@@ -363,7 +364,7 @@
sintheta = math.sin(theta_i)
while r_i < r0 + dr_i:
if rmax <= r_i:
- r_i += 0.5*dx
+ r_i += 0.5*dx
continue
y = r_i * costheta
x = r_i * sintheta
@@ -374,7 +375,7 @@
if img[pi, pj] != img[pi, pj]:
img[pi, pj] = 0.0
img[pi, pj] = field[i]
- r_i += 0.5*dx
+ r_i += 0.5*dx
theta_i += dthetamin
return img
@@ -412,7 +413,7 @@
cdef np.float64_t s2 = math.sqrt(2.0)
cdef np.float64_t xmax, ymax, xmin, ymin
nf = field.shape[0]
-
+
if input_img is None:
img = np.zeros((buff_size[0], buff_size[1]))
img[:] = np.nan
@@ -551,7 +552,7 @@
np.ndarray[np.int64_t, ndim=2] conn,
buff_size,
np.ndarray[np.float64_t, ndim=2] field,
- extents,
+ extents,
int index_offset = 0):
cdef np.ndarray[np.float64_t, ndim=3] img
img = np.zeros(buff_size, dtype="float64")
@@ -594,6 +595,8 @@
sampler = P1Sampler2D()
elif ndim == 2 and nvertices == 4:
sampler = Q1Sampler2D()
+ elif ndim == 2 and nvertices == 6:
+ sampler = T2Sampler2D()
else:
raise YTElementTypeNotRecognized(ndim, nvertices)
@@ -602,7 +605,7 @@
if buff_size[2] != 1:
raise RuntimeError("Slices of 2D datasets must be "
"perpendicular to the 'z' direction.")
-
+
# allocate temporary storage
vertices = <np.float64_t *> malloc(ndim * sizeof(np.float64_t) * nvertices)
field_vals = <np.float64_t *> malloc(sizeof(np.float64_t) * num_field_vals)
https://bitbucket.org/yt_analysis/yt/commits/fcdf0bcc26e5/
Changeset: fcdf0bcc26e5
Branch: yt
User: al007
Date: 2016-09-15 18:06:54+00:00
Summary: Update element_mapping description for second-order tris.
Affected #: 1 file
diff -r f9843e2a1748165775d26c9e458f88ff146a4424 -r fcdf0bcc26e582618ada707d021349691f5b7c21 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -830,8 +830,8 @@
@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.
+ # for canonical tris, we check whether the mapped_coords are between
+ # 0 and 1.
cdef int i
for i in range(2):
if (mapped_coord[i] < -self.inclusion_tol or
https://bitbucket.org/yt_analysis/yt/commits/e3698e365fb0/
Changeset: e3698e365fb0
Branch: yt
User: al007
Date: 2016-09-19 21:11:04+00:00
Summary: Unexpected EOF while parsing shape functions.
Affected #: 3 files
diff -r fcdf0bcc26e582618ada707d021349691f5b7c21 -r e3698e365fb093ab869ec17e9716a3ceade93ef9 yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -25,30 +25,3 @@
double* phys_x) nogil
-cdef void T2Function2D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil
-
-
-cdef void T2Jacobian2D(double* rcol,
- double* scol,
- double* x,
- double* vertices,
- double* phys_x) nogil
-
-
-cdef void W1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil
-
-
-cdef void W1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil
-
-
diff -r fcdf0bcc26e582618ada707d021349691f5b7c21 -r e3698e365fb093ab869ec17e9716a3ceade93ef9 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -65,60 +65,3 @@
scol[1] = -0.25*(1 - x[0])*vertices[1] + 0.25*(1 - x[0])*vertices[7] - 0.25*(1 + x[0])*vertices[3] + 0.25*(1 + x[0])*vertices[5];
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void T2Function2D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- fx[0] = (-x[0] + 2*pow(x[0], 2))*vertices[2] + (-x[1] + 2*pow(x[1], 2))*vertices[4] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[10] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[6] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[0] - phys_x[0] + 4*x[0]*x[1]*vertices[8];
- fx[1] = (-x[0] + 2*pow(x[0], 2))*vertices[3] + (-x[1] + 2*pow(x[1], 2))*vertices[5] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[11] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[7] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[1] - phys_x[1] + 4*x[0]*x[1]*vertices[9];
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void T2Jacobian2D(double* rcol,
- double* scol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- rcol[0] = (-1 + 4*x[0])*vertices[2] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 8*x[0] - 4*x[1])*vertices[6] + 4*x[1]*vertices[8] - 4*x[1]*vertices[10];
- scol[0] = (-1 + 4*x[1])*vertices[4] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 4*x[0] - 8*x[1])*vertices[10] - 4*x[0]*vertices[6] + 4*x[0]*vertices[8];
- rcol[1] = (-1 + 4*x[0])*vertices[3] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 8*x[0] - 4*x[1])*vertices[7] + 4*x[1]*vertices[9] - 4*x[1]*vertices[11];
- scol[1] = (-1 + 4*x[1])*vertices[5] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 4*x[0] - 8*x[1])*vertices[11] - 4*x[0]*vertices[7] + 4*x[0]*vertices[9];
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void W1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- fx[0] = 0.5*(1 - x[0] - x[1])*(1 - x[2])*vertices[0] + 0.5*(1 - x[0] - x[1])*(1 + x[2])*vertices[9] - phys_x[0] + 0.5*x[0]*(1 - x[2])*vertices[3] + 0.5*x[0]*(1 + x[2])*vertices[12] + 0.5*x[1]*(1 - x[2])*vertices[6] + 0.5*x[1]*(1 + x[2])*vertices[15];
- fx[1] = 0.5*(1 - x[0] - x[1])*(1 - x[2])*vertices[1] + 0.5*(1 - x[0] - x[1])*(1 + x[2])*vertices[10] - phys_x[1] + 0.5*x[0]*(1 - x[2])*vertices[4] + 0.5*x[0]*(1 + x[2])*vertices[13] + 0.5*x[1]*(1 - x[2])*vertices[7] + 0.5*x[1]*(1 + x[2])*vertices[16];
- fx[2] = 0.5*(1 - x[0] - x[1])*(1 - x[2])*vertices[2] + 0.5*(1 - x[0] - x[1])*(1 + x[2])*vertices[11] - phys_x[2] + 0.5*x[0]*(1 - x[2])*vertices[5] + 0.5*x[0]*(1 + x[2])*vertices[14] + 0.5*x[1]*(1 - x[2])*vertices[8] + 0.5*x[1]*(1 + x[2])*vertices[17];
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void W1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- rcol[0] = -0.5*(1 - x[2])*vertices[0] + 0.5*(1 - x[2])*vertices[3] - 0.5*(1 + x[2])*vertices[9] + 0.5*(1 + x[2])*vertices[12];
- scol[0] = -0.5*(1 - x[2])*vertices[0] + 0.5*(1 - x[2])*vertices[6] - 0.5*(1 + x[2])*vertices[9] + 0.5*(1 + x[2])*vertices[15];
- tcol[0] = -0.5*(1 - x[0] - x[1])*vertices[0] + 0.5*(1 - x[0] - x[1])*vertices[9] - 0.5*x[0]*vertices[3] + 0.5*x[0]*vertices[12] - 0.5*x[1]*vertices[6] + 0.5*x[1]*vertices[15];
- rcol[1] = -0.5*(1 - x[2])*vertices[1] + 0.5*(1 - x[2])*vertices[4] - 0.5*(1 + x[2])*vertices[10] + 0.5*(1 + x[2])*vertices[13];
- scol[1] = -0.5*(1 - x[2])*vertices[1] + 0.5*(1 - x[2])*vertices[7] - 0.5*(1 + x[2])*vertices[10] + 0.5*(1 + x[2])*vertices[16];
- tcol[1] = -0.5*(1 - x[0] - x[1])*vertices[1] + 0.5*(1 - x[0] - x[1])*vertices[10] - 0.5*x[0]*vertices[4] + 0.5*x[0]*vertices[13] - 0.5*x[1]*vertices[7] + 0.5*x[1]*vertices[16];
- rcol[2] = -0.5*(1 - x[2])*vertices[2] + 0.5*(1 - x[2])*vertices[5] - 0.5*(1 + x[2])*vertices[11] + 0.5*(1 + x[2])*vertices[14];
- scol[2] = -0.5*(1 - x[2])*vertices[2] + 0.5*(1 - x[2])*vertices[8] - 0.5*(1 + x[2])*vertices[11] + 0.5*(1 + x[2])*vertices[17];
- tcol[2] = -0.5*(1 - x[0] - x[1])*vertices[2] + 0.5*(1 - x[0] - x[1])*vertices[11] - 0.5*x[0]*vertices[5] + 0.5*x[0]*vertices[14] - 0.5*x[1]*vertices[8] + 0.5*x[1]*vertices[17];
-
-
diff -r fcdf0bcc26e582618ada707d021349691f5b7c21 -r e3698e365fb093ab869ec17e9716a3ceade93ef9 yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -49,3 +49,22 @@
4 * x[0] - 4 * x[0]**2 - 4 * x[0] * x[1],
4 * x[0] * x[1],
4 * x[1] - 4 * x[1]**2 - 4 * x[0] * x[1]]
+
+Tet10:
+ mesh_type: Tet2
+ num_dim: 3
+ num_vertices: 10
+ num_mapped_coords: 3
+ shape_functions: |
+ [1 - 3 * x[0] + 2 * x[0]**2 - 3 * x[1] + 2 * x[1]**2 - 3 * x[2] + 2 * x[2]**2 + 4 * x[0] * x[1] + 4 * x[0] * x[2] + 4 * x[1] * x[2],
+ -x[0] + 2 * x[0]**2,
+ -x[1] + 2 * x[1]**2,
+ -x[2] + 2 * x[2]**2,
+ 4 * x[0] - 4 * x[0]**2 - 4 * x[0] * x[1] - 4 * x[0] * x[2],
+ 4 * x[0] * x[1],
+ 4 * x[1] - 4 * x[1]**2 - 4 * x[1] * x[0] - 4 * x[1] * x[2],
+ 4 * x[2] - 4 * x[2]**2 - 4 * x[2] * x[0] - 4 * x[2] * x[0],
+ 4 * x[0] * x[2],
+ 4 * x[1] * x[2]]
+
+
\ No newline at end of file
https://bitbucket.org/yt_analysis/yt/commits/5ff4f827dc30/
Changeset: 5ff4f827dc30
Branch: yt
User: al007
Date: 2016-09-19 22:36:28+00:00
Summary: Implement visualization support for 2nd order tets.
Affected #: 7 files
diff -r e3698e365fb093ab869ec17e9716a3ceade93ef9 -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -25,3 +25,44 @@
double* phys_x) nogil
+cdef void Tet2Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Tet2Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void T2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void T2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void W1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void W1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
diff -r e3698e365fb093ab869ec17e9716a3ceade93ef9 -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -65,3 +65,92 @@
scol[1] = -0.25*(1 - x[0])*vertices[1] + 0.25*(1 - x[0])*vertices[7] - 0.25*(1 + x[0])*vertices[3] + 0.25*(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,
+ double* phys_x) nogil:
+ fx[0] = (-x[0] + 2*pow(x[0], 2))*vertices[3] + (-x[1] + 2*pow(x[1], 2))*vertices[6] + (-x[2] + 2*pow(x[2], 2))*vertices[9] + (4*x[0] - 4*x[0]*x[1] - 4*x[0]*x[2] - 4*pow(x[0], 2))*vertices[12] + (4*x[1] - 4*x[1]*x[0] - 4*x[1]*x[2] - 4*pow(x[1], 2))*vertices[18] + (4*x[2] - 4*x[2]*x[0] - 4*x[2]*x[1] - 4*pow(x[2], 2))*vertices[21] + (1 - 3*x[0] + 4*x[0]*x[1] + 4*x[0]*x[2] + 2*pow(x[0], 2) - 3*x[1] + 4*x[1]*x[2] + 2*pow(x[1], 2) - 3*x[2] + 2*pow(x[2], 2))*vertices[0] - phys_x[0] + 4*x[0]*x[1]*vertices[15] + 4*x[0]*x[2]*vertices[24] + 4*x[1]*x[2]*vertices[27];
+ fx[1] = (-x[0] + 2*pow(x[0], 2))*vertices[4] + (-x[1] + 2*pow(x[1], 2))*vertices[7] + (-x[2] + 2*pow(x[2], 2))*vertices[10] + (4*x[0] - 4*x[0]*x[1] - 4*x[0]*x[2] - 4*pow(x[0], 2))*vertices[13] + (4*x[1] - 4*x[1]*x[0] - 4*x[1]*x[2] - 4*pow(x[1], 2))*vertices[19] + (4*x[2] - 4*x[2]*x[0] - 4*x[2]*x[1] - 4*pow(x[2], 2))*vertices[22] + (1 - 3*x[0] + 4*x[0]*x[1] + 4*x[0]*x[2] + 2*pow(x[0], 2) - 3*x[1] + 4*x[1]*x[2] + 2*pow(x[1], 2) - 3*x[2] + 2*pow(x[2], 2))*vertices[1] - phys_x[1] + 4*x[0]*x[1]*vertices[16] + 4*x[0]*x[2]*vertices[25] + 4*x[1]*x[2]*vertices[28];
+ fx[2] = (-x[0] + 2*pow(x[0], 2))*vertices[5] + (-x[1] + 2*pow(x[1], 2))*vertices[8] + (-x[2] + 2*pow(x[2], 2))*vertices[11] + (4*x[0] - 4*x[0]*x[1] - 4*x[0]*x[2] - 4*pow(x[0], 2))*vertices[14] + (4*x[1] - 4*x[1]*x[0] - 4*x[1]*x[2] - 4*pow(x[1], 2))*vertices[20] + (4*x[2] - 4*x[2]*x[0] - 4*x[2]*x[1] - 4*pow(x[2], 2))*vertices[23] + (1 - 3*x[0] + 4*x[0]*x[1] + 4*x[0]*x[2] + 2*pow(x[0], 2) - 3*x[1] + 4*x[1]*x[2] + 2*pow(x[1], 2) - 3*x[2] + 2*pow(x[2], 2))*vertices[2] - phys_x[2] + 4*x[0]*x[1]*vertices[17] + 4*x[0]*x[2]*vertices[26] + 4*x[1]*x[2]*vertices[29];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Tet2Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = (-1 + 4*x[0])*vertices[3] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[0] + (4 - 8*x[0] - 4*x[1] - 4*x[2])*vertices[12] + 4*x[1]*vertices[15] - 4*x[1]*vertices[18] - 4*x[2]*vertices[21] + 4*x[2]*vertices[24];
+ scol[0] = (-1 + 4*x[1])*vertices[6] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[0] + (4 - 4*x[0] - 8*x[1] - 4*x[2])*vertices[18] - 4*x[0]*vertices[12] + 4*x[0]*vertices[15] - 4*x[2]*vertices[21] + 4*x[2]*vertices[27];
+ tcol[0] = (-1 + 4*x[2])*vertices[9] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[0] + (4 - 4*x[0] - 4*x[1] - 8*x[2])*vertices[21] - 4*x[0]*vertices[12] + 4*x[0]*vertices[24] - 4*x[1]*vertices[18] + 4*x[1]*vertices[27];
+ rcol[1] = (-1 + 4*x[0])*vertices[4] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[1] + (4 - 8*x[0] - 4*x[1] - 4*x[2])*vertices[13] + 4*x[1]*vertices[16] - 4*x[1]*vertices[19] - 4*x[2]*vertices[22] + 4*x[2]*vertices[25];
+ scol[1] = (-1 + 4*x[1])*vertices[7] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[1] + (4 - 4*x[0] - 8*x[1] - 4*x[2])*vertices[19] - 4*x[0]*vertices[13] + 4*x[0]*vertices[16] - 4*x[2]*vertices[22] + 4*x[2]*vertices[28];
+ tcol[1] = (-1 + 4*x[2])*vertices[10] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[1] + (4 - 4*x[0] - 4*x[1] - 8*x[2])*vertices[22] - 4*x[0]*vertices[13] + 4*x[0]*vertices[25] - 4*x[1]*vertices[19] + 4*x[1]*vertices[28];
+ rcol[2] = (-1 + 4*x[0])*vertices[5] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[2] + (4 - 8*x[0] - 4*x[1] - 4*x[2])*vertices[14] + 4*x[1]*vertices[17] - 4*x[1]*vertices[20] - 4*x[2]*vertices[23] + 4*x[2]*vertices[26];
+ scol[2] = (-1 + 4*x[1])*vertices[8] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[2] + (4 - 4*x[0] - 8*x[1] - 4*x[2])*vertices[20] - 4*x[0]*vertices[14] + 4*x[0]*vertices[17] - 4*x[2]*vertices[23] + 4*x[2]*vertices[29];
+ tcol[2] = (-1 + 4*x[2])*vertices[11] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[2] + (4 - 4*x[0] - 4*x[1] - 8*x[2])*vertices[23] - 4*x[0]*vertices[14] + 4*x[0]*vertices[26] - 4*x[1]*vertices[20] + 4*x[1]*vertices[29];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void T2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = (-x[0] + 2*pow(x[0], 2))*vertices[2] + (-x[1] + 2*pow(x[1], 2))*vertices[4] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[10] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[6] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[0] - phys_x[0] + 4*x[0]*x[1]*vertices[8];
+ fx[1] = (-x[0] + 2*pow(x[0], 2))*vertices[3] + (-x[1] + 2*pow(x[1], 2))*vertices[5] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[11] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[7] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[1] - phys_x[1] + 4*x[0]*x[1]*vertices[9];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void T2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = (-1 + 4*x[0])*vertices[2] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 8*x[0] - 4*x[1])*vertices[6] + 4*x[1]*vertices[8] - 4*x[1]*vertices[10];
+ scol[0] = (-1 + 4*x[1])*vertices[4] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 4*x[0] - 8*x[1])*vertices[10] - 4*x[0]*vertices[6] + 4*x[0]*vertices[8];
+ rcol[1] = (-1 + 4*x[0])*vertices[3] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 8*x[0] - 4*x[1])*vertices[7] + 4*x[1]*vertices[9] - 4*x[1]*vertices[11];
+ scol[1] = (-1 + 4*x[1])*vertices[5] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 4*x[0] - 8*x[1])*vertices[11] - 4*x[0]*vertices[7] + 4*x[0]*vertices[9];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void W1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = 0.5*(1 - x[0] - x[1])*(1 - x[2])*vertices[0] + 0.5*(1 - x[0] - x[1])*(1 + x[2])*vertices[9] - phys_x[0] + 0.5*x[0]*(1 - x[2])*vertices[3] + 0.5*x[0]*(1 + x[2])*vertices[12] + 0.5*x[1]*(1 - x[2])*vertices[6] + 0.5*x[1]*(1 + x[2])*vertices[15];
+ fx[1] = 0.5*(1 - x[0] - x[1])*(1 - x[2])*vertices[1] + 0.5*(1 - x[0] - x[1])*(1 + x[2])*vertices[10] - phys_x[1] + 0.5*x[0]*(1 - x[2])*vertices[4] + 0.5*x[0]*(1 + x[2])*vertices[13] + 0.5*x[1]*(1 - x[2])*vertices[7] + 0.5*x[1]*(1 + x[2])*vertices[16];
+ fx[2] = 0.5*(1 - x[0] - x[1])*(1 - x[2])*vertices[2] + 0.5*(1 - x[0] - x[1])*(1 + x[2])*vertices[11] - phys_x[2] + 0.5*x[0]*(1 - x[2])*vertices[5] + 0.5*x[0]*(1 + x[2])*vertices[14] + 0.5*x[1]*(1 - x[2])*vertices[8] + 0.5*x[1]*(1 + x[2])*vertices[17];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void W1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.5*(1 - x[2])*vertices[0] + 0.5*(1 - x[2])*vertices[3] - 0.5*(1 + x[2])*vertices[9] + 0.5*(1 + x[2])*vertices[12];
+ scol[0] = -0.5*(1 - x[2])*vertices[0] + 0.5*(1 - x[2])*vertices[6] - 0.5*(1 + x[2])*vertices[9] + 0.5*(1 + x[2])*vertices[15];
+ tcol[0] = -0.5*(1 - x[0] - x[1])*vertices[0] + 0.5*(1 - x[0] - x[1])*vertices[9] - 0.5*x[0]*vertices[3] + 0.5*x[0]*vertices[12] - 0.5*x[1]*vertices[6] + 0.5*x[1]*vertices[15];
+ rcol[1] = -0.5*(1 - x[2])*vertices[1] + 0.5*(1 - x[2])*vertices[4] - 0.5*(1 + x[2])*vertices[10] + 0.5*(1 + x[2])*vertices[13];
+ scol[1] = -0.5*(1 - x[2])*vertices[1] + 0.5*(1 - x[2])*vertices[7] - 0.5*(1 + x[2])*vertices[10] + 0.5*(1 + x[2])*vertices[16];
+ tcol[1] = -0.5*(1 - x[0] - x[1])*vertices[1] + 0.5*(1 - x[0] - x[1])*vertices[10] - 0.5*x[0]*vertices[4] + 0.5*x[0]*vertices[13] - 0.5*x[1]*vertices[7] + 0.5*x[1]*vertices[16];
+ rcol[2] = -0.5*(1 - x[2])*vertices[2] + 0.5*(1 - x[2])*vertices[5] - 0.5*(1 + x[2])*vertices[11] + 0.5*(1 + x[2])*vertices[14];
+ scol[2] = -0.5*(1 - x[2])*vertices[2] + 0.5*(1 - x[2])*vertices[8] - 0.5*(1 + x[2])*vertices[11] + 0.5*(1 + x[2])*vertices[17];
+ tcol[2] = -0.5*(1 - x[0] - x[1])*vertices[2] + 0.5*(1 - x[0] - x[1])*vertices[11] - 0.5*x[0]*vertices[5] + 0.5*x[0]*vertices[14] - 0.5*x[1]*vertices[8] + 0.5*x[1]*vertices[17];
+
+
diff -r e3698e365fb093ab869ec17e9716a3ceade93ef9 -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -235,3 +235,17 @@
double* vals) nogil
cdef int check_inside(self, double* mapped_coord) nogil
+
+cdef class Tet2Sampler3D(NonlinearSolveSampler3D):
+
+ 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
diff -r e3698e365fb093ab869ec17e9716a3ceade93ef9 -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -27,7 +27,9 @@
W1Function3D, \
W1Jacobian3D, \
T2Function2D, \
- T2Jacobian2D
+ T2Jacobian2D, \
+ Tet2Function3D, \
+ Tet2Jacobian3D
cdef extern from "platform_dep.h":
double fmax(double x, double y) nogil
@@ -839,6 +841,65 @@
return 0
return 1
+cdef class Tet2Sampler3D(NonlinearSolveSampler3D):
+
+ '''
+
+ This implements sampling inside a 3D, quadratic, tetrahedral mesh
+ element. Note that this implementation uses canonical coordinates.
+
+ '''
+
+ def __init__(self):
+ super(Tet2Sampler3D, self).__init__()
+ self.num_mapped_coords = 3
+ self.dim = 3
+ self.func = Tet2Function3D
+ self.jac = Tet2Jacobian3D
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
+ cdef double[10] phi
+ cdef int i
+ cdef double return_value = 0
+
+ phi[0] = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
+ 2 * coord[1]**2 - 3 * coord[2] + 2 * coord[2]**2 + \
+ 4 * coord[0] * coord[1] + 4 * coord[0] * coord[2] + \
+ 4 * coord[1] * coord[2]
+ phi[1] = -coord[0] + 2 * coord[0]**2
+ phi[2] = -coord[1] + 2 * coord[1]**2
+ phi[3] = -coord[2] + 2 * coord[2]**2
+ phi[4] = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1] - \
+ 4 * coord[0] * coord[2]
+ phi[5] = 4 * coord[0] * coord[1]
+ phi[6] = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1] - \
+ 4 * coord[1] * coord[2]
+ phi[7] = 4 * coord[2] - 4 * coord[2]**2 - 4 * coord[2] * coord[0] - \
+ 4 * coord[2] * coord[1]
+ phi[8] = 4 * coord[0] * coord[2]
+ phi[9] = 4 * coord[1] * coord[2]
+
+ for i in range(10):
+ return_value += phi[i] * vals[i]
+
+ return return_value
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_inside(self, double* mapped_coord) nogil:
+ # for canonical tets, we check whether the mapped_coords are between
+ # 0 and 1.
+ cdef int i
+ for i in range(3):
+ if (mapped_coord[i] < -self.inclusion_tol or
+ mapped_coord[i] - 1.0 > self.inclusion_tol):
+ return 0
+ return 1
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -959,3 +1020,20 @@
<double*> physical_x.data)
return val
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+def test_tet2_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 = Tet2Sampler3D()
+
+ val = sampler.sample_at_real_point(<double*> vertices.data,
+ <double*> field_values.data,
+ <double*> physical_x.data)
+
+ return val
diff -r e3698e365fb093ab869ec17e9716a3ceade93ef9 -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -31,7 +31,8 @@
P1Sampler2D, \
Q1Sampler2D, \
W1Sampler3D, \
- T2Sampler2D
+ T2Sampler2D, \
+ Tet2Sampler3D
cdef extern from "pixelization_constants.h":
enum:
@@ -597,6 +598,8 @@
sampler = Q1Sampler2D()
elif ndim == 2 and nvertices == 6:
sampler = T2Sampler2D()
+ elif ndim == 3 and nvertices == 10:
+ sampler = Tet2Sampler3D()
else:
raise YTElementTypeNotRecognized(ndim, nvertices)
diff -r e3698e365fb093ab869ec17e9716a3ceade93ef9 -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 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
@@ -23,7 +23,8 @@
test_quad_sampler, \
test_hex20_sampler, \
test_wedge_sampler, \
- test_tri2_sampler
+ test_tri2_sampler, \
+ test_tet2_sampler
def check_all_vertices(sampler, vertices, field_values):
@@ -33,6 +34,8 @@
for i in range(NV):
x = vertices[i]
val = sampler(vertices, field_values, x)
+ print("Actual equals " + str(val))
+ print("Desired equals " + str(field_values[i]))
assert_almost_equal(val, field_values[i])
@@ -141,3 +144,21 @@
field_values = np.array([15., 37., 49., 32., 46., 24.])
check_all_vertices(test_tri2_sampler, vertices, field_values)
+
+
+def test_Tet2Sampler3D():
+
+ vertices = np.array([[0.3 , -0.4 , 0.6] ,
+ [1.7 , -0.7 , 0.8] ,
+ [0.4 , 1.2 , 0.4] ,
+ [0.4 , -0.2 , 2.0] ,
+ [1.0 , -0.55 , 0.7] ,
+ [1.05 , 0.25 , 0.6] ,
+ [0.35 , 0.4 , 0.5] ,
+ [0.35 , -0.3 , 1.3] ,
+ [1.05 , -0.45 , 1.4] ,
+ [0.4 , 0.5 , 1.2]])
+
+ field_values = np.array([15., 37., 49., 24., 30., 44., 20., 17., 32., 36.])
+
+ check_all_vertices(test_tet2_sampler, vertices, field_values)
diff -r e3698e365fb093ab869ec17e9716a3ceade93ef9 -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -63,8 +63,6 @@
4 * x[0] - 4 * x[0]**2 - 4 * x[0] * x[1] - 4 * x[0] * x[2],
4 * x[0] * x[1],
4 * x[1] - 4 * x[1]**2 - 4 * x[1] * x[0] - 4 * x[1] * x[2],
- 4 * x[2] - 4 * x[2]**2 - 4 * x[2] * x[0] - 4 * x[2] * x[0],
+ 4 * x[2] - 4 * x[2]**2 - 4 * x[2] * x[0] - 4 * x[2] * x[1],
4 * x[0] * x[2],
4 * x[1] * x[2]]
-
-
\ No newline at end of file
https://bitbucket.org/yt_analysis/yt/commits/312b7acc1917/
Changeset: 312b7acc1917
Branch: yt
User: al007
Date: 2016-09-20 22:14:17+00:00
Summary: Add answer test for second order tris.
Affected #: 3 files
diff -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 -r 312b7acc191727206d6f5969c6b11ad87460f233 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -35,10 +35,10 @@
- yt/frontends/owls_subfind/tests/test_outputs.py
- yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g5
- yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g42
-
+
local_owls_000:
- yt/frontends/owls/tests/test_outputs.py
-
+
local_pw_006:
- yt/visualization/tests/test_plotwindow.py:test_attributes
- yt/visualization/tests/test_plotwindow.py:test_attributes_wt
@@ -46,10 +46,10 @@
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_answers
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_filter
- yt/visualization/tests/test_particle_plot.py:test_particle_phase_answers
-
+
local_tipsy_001:
- yt/frontends/tipsy/tests/test_outputs.py
-
+
local_varia_003:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
@@ -57,13 +57,14 @@
- yt/analysis_modules/photon_simulator/tests/test_sloshing.py
- yt/visualization/volume_rendering/tests/test_vr_orientation.py
- yt/visualization/volume_rendering/tests/test_mesh_render.py
+ - yt/visualization/test_mesh_slices.py
local_orion_000:
- yt/frontends/boxlib/tests/test_orion.py
-
+
local_ramses_000:
- yt/frontends/ramses/tests/test_outputs.py
-
+
local_ytdata_000:
- yt/frontends/ytdata
diff -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 -r 312b7acc191727206d6f5969c6b11ad87460f233 yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -714,7 +714,7 @@
def compare(self, new_result, old_result):
compare_image_lists(new_result, old_result, self.decimals)
-
+
class PlotWindowAttributeTest(AnswerTestingTest):
_type_name = "PlotWindowAttribute"
_attrs = ('plot_type', 'plot_field', 'plot_axis', 'attr_name', 'attr_args',
@@ -758,7 +758,7 @@
_type_name = "PhasePlotAttribute"
_attrs = ('plot_type', 'x_field', 'y_field', 'z_field',
'attr_name', 'attr_args')
- def __init__(self, ds_fn, x_field, y_field, z_field,
+ def __init__(self, ds_fn, x_field, y_field, z_field,
attr_name, attr_args, decimals, plot_type='PhasePlot'):
super(PhasePlotAttributeTest, self).__init__(ds_fn)
self.data_source = self.ds.all_data()
@@ -771,7 +771,7 @@
self.attr_args = attr_args
self.decimals = decimals
- def create_plot(self, data_source, x_field, y_field, z_field,
+ def create_plot(self, data_source, x_field, y_field, z_field,
plot_type, plot_kwargs=None):
# plot_type should be a string
# plot_kwargs should be a dict
@@ -928,7 +928,7 @@
return ftrue
else:
return ffalse
-
+
def requires_ds(ds_fn, big_data = False, file_check = False):
def ffalse(func):
return lambda: None
diff -r 5ff4f827dc3049393cd653cfd20e0ca648c23d43 -r 312b7acc191727206d6f5969c6b11ad87460f233 yt/visualization/tests/test_mesh_slices.py
--- a/yt/visualization/tests/test_mesh_slices.py
+++ b/yt/visualization/tests/test_mesh_slices.py
@@ -18,6 +18,10 @@
import yt
from yt.testing import fake_tetrahedral_ds
from yt.testing import fake_hexahedral_ds
+from yt.utilities.answer_testing.framework import \
+ requires_ds, \
+ data_dir_load, \
+ GenericImageTest
def setup():
@@ -25,6 +29,23 @@
from yt.config import ytcfg
ytcfg["yt", "__withintesting"] = "True"
+def compare(ds, field, test_prefix, decimals=12):
+ def slice_image(filename_prefix):
+ sl = yt.SlicePlot(ds, 'z', field)
+ return sl.save(filename_prefix)
+
+ slice_image.__name__ = "slice_{}".format(test_prefix)
+ test = GenericImageTest(ds, slice_image, decimals)
+ test.prefix = test_prefix
+ return test
+
+tri2 = "MOOSE_sample_data/RZ_p_no_parts_do_nothing_bcs_cone_out.e"
+
+ at requires_ds(tri2)
+def test_tri2():
+ ds = data_dir_load(tri2)
+ for field in ds.field_list:
+ yield compare(ds, field, "answers_tri2_%s_%s" % (field[0], field[1]))
def test_mesh_slices():
# Perform I/O in safe place instead of yt main dir
https://bitbucket.org/yt_analysis/yt/commits/174177f6c363/
Changeset: 174177f6c363
Branch: yt
User: al007
Date: 2016-09-21 18:40:25+00:00
Summary: Load correct time step for answer test.
Affected #: 1 file
diff -r 312b7acc191727206d6f5969c6b11ad87460f233 -r 174177f6c36314e498f9a6f4f8d4fc87bb0e4ba7 yt/visualization/tests/test_mesh_slices.py
--- a/yt/visualization/tests/test_mesh_slices.py
+++ b/yt/visualization/tests/test_mesh_slices.py
@@ -32,7 +32,9 @@
def compare(ds, field, test_prefix, decimals=12):
def slice_image(filename_prefix):
sl = yt.SlicePlot(ds, 'z', field)
- return sl.save(filename_prefix)
+ sl.set_log('all', False)
+ image_file = sl.save(filename_prefix)
+ return image_file
slice_image.__name__ = "slice_{}".format(test_prefix)
test = GenericImageTest(ds, slice_image, decimals)
@@ -43,7 +45,7 @@
@requires_ds(tri2)
def test_tri2():
- ds = data_dir_load(tri2)
+ ds = data_dir_load(tri2, kwargs={'step':-1})
for field in ds.field_list:
yield compare(ds, field, "answers_tri2_%s_%s" % (field[0], field[1]))
https://bitbucket.org/yt_analysis/yt/commits/ebc0cd2da2a0/
Changeset: ebc0cd2da2a0
Branch: yt
User: al007
Date: 2016-09-21 18:52:01+00:00
Summary: Fix answer test location in tests.yaml and remove print statements.
Affected #: 2 files
diff -r 174177f6c36314e498f9a6f4f8d4fc87bb0e4ba7 -r ebc0cd2da2a0b55a275254f371872b275e1a3d82 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -57,7 +57,7 @@
- yt/analysis_modules/photon_simulator/tests/test_sloshing.py
- yt/visualization/volume_rendering/tests/test_vr_orientation.py
- yt/visualization/volume_rendering/tests/test_mesh_render.py
- - yt/visualization/test_mesh_slices.py
+ - yt/visualization/tests/test_mesh_slices.py
local_orion_000:
- yt/frontends/boxlib/tests/test_orion.py
diff -r 174177f6c36314e498f9a6f4f8d4fc87bb0e4ba7 -r ebc0cd2da2a0b55a275254f371872b275e1a3d82 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
@@ -34,8 +34,6 @@
for i in range(NV):
x = vertices[i]
val = sampler(vertices, field_values, x)
- print("Actual equals " + str(val))
- print("Desired equals " + str(field_values[i]))
assert_almost_equal(val, field_values[i])
https://bitbucket.org/yt_analysis/yt/commits/ef1826dd3249/
Changeset: ef1826dd3249
Branch: yt
User: al007
Date: 2016-09-21 20:29:48+00:00
Summary: Update data set location.
Affected #: 1 file
diff -r ebc0cd2da2a0b55a275254f371872b275e1a3d82 -r ef1826dd3249232d5874131f51292818493e71a7 yt/visualization/tests/test_mesh_slices.py
--- a/yt/visualization/tests/test_mesh_slices.py
+++ b/yt/visualization/tests/test_mesh_slices.py
@@ -41,7 +41,7 @@
test.prefix = test_prefix
return test
-tri2 = "MOOSE_sample_data/RZ_p_no_parts_do_nothing_bcs_cone_out.e"
+tri2 = "SecondOrderTris/RZ_p_no_parts_do_nothing_bcs_cone_out.e"
@requires_ds(tri2)
def test_tri2():
https://bitbucket.org/yt_analysis/yt/commits/0150db45e998/
Changeset: 0150db45e998
Branch: yt
User: al007
Date: 2016-09-22 02:21:22+00:00
Summary: Fix up tests.yaml.
Affected #: 1 file
diff -r ef1826dd3249232d5874131f51292818493e71a7 -r 0150db45e998c536eeb4bb95bd2aa2fff9ea23ba tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -50,14 +50,14 @@
local_tipsy_001:
- yt/frontends/tipsy/tests/test_outputs.py
- local_varia_003:
+ local_varia_004:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
- yt/analysis_modules/photon_simulator/tests/test_sloshing.py
- yt/visualization/volume_rendering/tests/test_vr_orientation.py
- yt/visualization/volume_rendering/tests/test_mesh_render.py
- - yt/visualization/tests/test_mesh_slices.py
+ - yt/visualization/tests/test_mesh_slices.py:test_tri2
local_orion_000:
- yt/frontends/boxlib/tests/test_orion.py
https://bitbucket.org/yt_analysis/yt/commits/b1e7a7e87d3f/
Changeset: b1e7a7e87d3f
Branch: yt
User: ngoldbaum
Date: 2016-09-22 03:32:22+00:00
Summary: Merged in al007/yt (pull request #2378)
Add support for second order triangles in unstructured meshes
Affected #: 13 files
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -35,10 +35,10 @@
- yt/frontends/owls_subfind/tests/test_outputs.py
- yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g5
- yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g42
-
+
local_owls_000:
- yt/frontends/owls/tests/test_outputs.py
-
+
local_pw_006:
- yt/visualization/tests/test_plotwindow.py:test_attributes
- yt/visualization/tests/test_plotwindow.py:test_attributes_wt
@@ -46,24 +46,25 @@
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_answers
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_filter
- yt/visualization/tests/test_particle_plot.py:test_particle_phase_answers
-
+
local_tipsy_001:
- yt/frontends/tipsy/tests/test_outputs.py
-
- local_varia_003:
+
+ local_varia_004:
- yt/analysis_modules/radmc3d_export
- yt/frontends/moab/tests/test_c5.py
- yt/analysis_modules/photon_simulator/tests/test_spectra.py
- yt/analysis_modules/photon_simulator/tests/test_sloshing.py
- yt/visualization/volume_rendering/tests/test_vr_orientation.py
- yt/visualization/volume_rendering/tests/test_mesh_render.py
+ - yt/visualization/tests/test_mesh_slices.py:test_tri2
local_orion_000:
- yt/frontends/boxlib/tests/test_orion.py
-
+
local_ramses_000:
- yt/frontends/ramses/tests/test_outputs.py
-
+
local_ytdata_000:
- yt/frontends/ytdata
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -172,4 +172,3 @@
@property
def period(self):
return self.ds.domain_width
-
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/geometry/geometry_handler.py
--- a/yt/geometry/geometry_handler.py
+++ b/yt/geometry/geometry_handler.py
@@ -270,7 +270,7 @@
else:
tr = func(self)
if self._cache:
-
+
setattr(self, n, tr)
return tr
return property(cached_func)
@@ -397,7 +397,9 @@
@cached_property
def fcoords_vertex(self):
- ci = np.empty((self.data_size, 8, 3), dtype='float64')
+ nodes_per_elem = self.dobj.index.meshes[0].connectivity_indices.shape[1]
+ dim = self.dobj.ds.dimensionality
+ ci = np.empty((self.data_size, nodes_per_elem, dim), dtype='float64')
ci = YTArray(ci, input_units = "code_length",
registry = self.dobj.ds.unit_registry)
if self.data_size == 0: return ci
@@ -426,10 +428,10 @@
def __iter__(self):
return self
-
+
def __next__(self):
return self.next()
-
+
def next(self):
if len(self.queue) == 0:
for i in range(self.max_length):
@@ -445,4 +447,3 @@
g = self.queue.pop(0)
g._initialize_cache(self.cache.pop(g.id, {}))
return g
-
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/utilities/answer_testing/framework.py
--- a/yt/utilities/answer_testing/framework.py
+++ b/yt/utilities/answer_testing/framework.py
@@ -714,7 +714,7 @@
def compare(self, new_result, old_result):
compare_image_lists(new_result, old_result, self.decimals)
-
+
class PlotWindowAttributeTest(AnswerTestingTest):
_type_name = "PlotWindowAttribute"
_attrs = ('plot_type', 'plot_field', 'plot_axis', 'attr_name', 'attr_args',
@@ -758,7 +758,7 @@
_type_name = "PhasePlotAttribute"
_attrs = ('plot_type', 'x_field', 'y_field', 'z_field',
'attr_name', 'attr_args')
- def __init__(self, ds_fn, x_field, y_field, z_field,
+ def __init__(self, ds_fn, x_field, y_field, z_field,
attr_name, attr_args, decimals, plot_type='PhasePlot'):
super(PhasePlotAttributeTest, self).__init__(ds_fn)
self.data_source = self.ds.all_data()
@@ -771,7 +771,7 @@
self.attr_args = attr_args
self.decimals = decimals
- def create_plot(self, data_source, x_field, y_field, z_field,
+ def create_plot(self, data_source, x_field, y_field, z_field,
plot_type, plot_kwargs=None):
# plot_type should be a string
# plot_kwargs should be a dict
@@ -928,7 +928,7 @@
return ftrue
else:
return ffalse
-
+
def requires_ds(ds_fn, big_data = False, file_check = False):
def ffalse(func):
return lambda: None
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac 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,33 @@
double* phys_x) nogil
+cdef void Tet2Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Tet2Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void T2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void T2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
cdef void W1Function3D(double* fx,
double* x,
double* vertices,
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -5,7 +5,7 @@
cimport cython
-
+from libc.math cimport pow
@cython.boundscheck(False)
@@ -68,6 +68,63 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+cdef void Tet2Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = (-x[0] + 2*pow(x[0], 2))*vertices[3] + (-x[1] + 2*pow(x[1], 2))*vertices[6] + (-x[2] + 2*pow(x[2], 2))*vertices[9] + (4*x[0] - 4*x[0]*x[1] - 4*x[0]*x[2] - 4*pow(x[0], 2))*vertices[12] + (4*x[1] - 4*x[1]*x[0] - 4*x[1]*x[2] - 4*pow(x[1], 2))*vertices[18] + (4*x[2] - 4*x[2]*x[0] - 4*x[2]*x[1] - 4*pow(x[2], 2))*vertices[21] + (1 - 3*x[0] + 4*x[0]*x[1] + 4*x[0]*x[2] + 2*pow(x[0], 2) - 3*x[1] + 4*x[1]*x[2] + 2*pow(x[1], 2) - 3*x[2] + 2*pow(x[2], 2))*vertices[0] - phys_x[0] + 4*x[0]*x[1]*vertices[15] + 4*x[0]*x[2]*vertices[24] + 4*x[1]*x[2]*vertices[27];
+ fx[1] = (-x[0] + 2*pow(x[0], 2))*vertices[4] + (-x[1] + 2*pow(x[1], 2))*vertices[7] + (-x[2] + 2*pow(x[2], 2))*vertices[10] + (4*x[0] - 4*x[0]*x[1] - 4*x[0]*x[2] - 4*pow(x[0], 2))*vertices[13] + (4*x[1] - 4*x[1]*x[0] - 4*x[1]*x[2] - 4*pow(x[1], 2))*vertices[19] + (4*x[2] - 4*x[2]*x[0] - 4*x[2]*x[1] - 4*pow(x[2], 2))*vertices[22] + (1 - 3*x[0] + 4*x[0]*x[1] + 4*x[0]*x[2] + 2*pow(x[0], 2) - 3*x[1] + 4*x[1]*x[2] + 2*pow(x[1], 2) - 3*x[2] + 2*pow(x[2], 2))*vertices[1] - phys_x[1] + 4*x[0]*x[1]*vertices[16] + 4*x[0]*x[2]*vertices[25] + 4*x[1]*x[2]*vertices[28];
+ fx[2] = (-x[0] + 2*pow(x[0], 2))*vertices[5] + (-x[1] + 2*pow(x[1], 2))*vertices[8] + (-x[2] + 2*pow(x[2], 2))*vertices[11] + (4*x[0] - 4*x[0]*x[1] - 4*x[0]*x[2] - 4*pow(x[0], 2))*vertices[14] + (4*x[1] - 4*x[1]*x[0] - 4*x[1]*x[2] - 4*pow(x[1], 2))*vertices[20] + (4*x[2] - 4*x[2]*x[0] - 4*x[2]*x[1] - 4*pow(x[2], 2))*vertices[23] + (1 - 3*x[0] + 4*x[0]*x[1] + 4*x[0]*x[2] + 2*pow(x[0], 2) - 3*x[1] + 4*x[1]*x[2] + 2*pow(x[1], 2) - 3*x[2] + 2*pow(x[2], 2))*vertices[2] - phys_x[2] + 4*x[0]*x[1]*vertices[17] + 4*x[0]*x[2]*vertices[26] + 4*x[1]*x[2]*vertices[29];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Tet2Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = (-1 + 4*x[0])*vertices[3] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[0] + (4 - 8*x[0] - 4*x[1] - 4*x[2])*vertices[12] + 4*x[1]*vertices[15] - 4*x[1]*vertices[18] - 4*x[2]*vertices[21] + 4*x[2]*vertices[24];
+ scol[0] = (-1 + 4*x[1])*vertices[6] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[0] + (4 - 4*x[0] - 8*x[1] - 4*x[2])*vertices[18] - 4*x[0]*vertices[12] + 4*x[0]*vertices[15] - 4*x[2]*vertices[21] + 4*x[2]*vertices[27];
+ tcol[0] = (-1 + 4*x[2])*vertices[9] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[0] + (4 - 4*x[0] - 4*x[1] - 8*x[2])*vertices[21] - 4*x[0]*vertices[12] + 4*x[0]*vertices[24] - 4*x[1]*vertices[18] + 4*x[1]*vertices[27];
+ rcol[1] = (-1 + 4*x[0])*vertices[4] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[1] + (4 - 8*x[0] - 4*x[1] - 4*x[2])*vertices[13] + 4*x[1]*vertices[16] - 4*x[1]*vertices[19] - 4*x[2]*vertices[22] + 4*x[2]*vertices[25];
+ scol[1] = (-1 + 4*x[1])*vertices[7] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[1] + (4 - 4*x[0] - 8*x[1] - 4*x[2])*vertices[19] - 4*x[0]*vertices[13] + 4*x[0]*vertices[16] - 4*x[2]*vertices[22] + 4*x[2]*vertices[28];
+ tcol[1] = (-1 + 4*x[2])*vertices[10] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[1] + (4 - 4*x[0] - 4*x[1] - 8*x[2])*vertices[22] - 4*x[0]*vertices[13] + 4*x[0]*vertices[25] - 4*x[1]*vertices[19] + 4*x[1]*vertices[28];
+ rcol[2] = (-1 + 4*x[0])*vertices[5] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[2] + (4 - 8*x[0] - 4*x[1] - 4*x[2])*vertices[14] + 4*x[1]*vertices[17] - 4*x[1]*vertices[20] - 4*x[2]*vertices[23] + 4*x[2]*vertices[26];
+ scol[2] = (-1 + 4*x[1])*vertices[8] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[2] + (4 - 4*x[0] - 8*x[1] - 4*x[2])*vertices[20] - 4*x[0]*vertices[14] + 4*x[0]*vertices[17] - 4*x[2]*vertices[23] + 4*x[2]*vertices[29];
+ tcol[2] = (-1 + 4*x[2])*vertices[11] + (-3 + 4*x[0] + 4*x[1] + 4*x[2])*vertices[2] + (4 - 4*x[0] - 4*x[1] - 8*x[2])*vertices[23] - 4*x[0]*vertices[14] + 4*x[0]*vertices[26] - 4*x[1]*vertices[20] + 4*x[1]*vertices[29];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void T2Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = (-x[0] + 2*pow(x[0], 2))*vertices[2] + (-x[1] + 2*pow(x[1], 2))*vertices[4] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[10] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[6] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[0] - phys_x[0] + 4*x[0]*x[1]*vertices[8];
+ fx[1] = (-x[0] + 2*pow(x[0], 2))*vertices[3] + (-x[1] + 2*pow(x[1], 2))*vertices[5] + (-4*x[0]*x[1] + 4*x[1] - 4*pow(x[1], 2))*vertices[11] + (4*x[0] - 4*x[0]*x[1] - 4*pow(x[0], 2))*vertices[7] + (1 - 3*x[0] + 4*x[0]*x[1] + 2*pow(x[0], 2) - 3*x[1] + 2*pow(x[1], 2))*vertices[1] - phys_x[1] + 4*x[0]*x[1]*vertices[9];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void T2Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = (-1 + 4*x[0])*vertices[2] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 8*x[0] - 4*x[1])*vertices[6] + 4*x[1]*vertices[8] - 4*x[1]*vertices[10];
+ scol[0] = (-1 + 4*x[1])*vertices[4] + (-3 + 4*x[0] + 4*x[1])*vertices[0] + (4 - 4*x[0] - 8*x[1])*vertices[10] - 4*x[0]*vertices[6] + 4*x[0]*vertices[8];
+ rcol[1] = (-1 + 4*x[0])*vertices[3] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 8*x[0] - 4*x[1])*vertices[7] + 4*x[1]*vertices[9] - 4*x[1]*vertices[11];
+ scol[1] = (-1 + 4*x[1])*vertices[5] + (-3 + 4*x[0] + 4*x[1])*vertices[1] + (4 - 4*x[0] - 8*x[1])*vertices[11] - 4*x[0]*vertices[7] + 4*x[0]*vertices[9];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
cdef void W1Function3D(double* fx,
double* x,
double* vertices,
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -13,7 +13,7 @@
cdef int num_mapped_coords
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -21,7 +21,7 @@
cdef double sample_at_unit_point(self,
double* coord,
double* vals) nogil
-
+
cdef double sample_at_real_point(self,
double* vertices,
@@ -36,7 +36,7 @@
cdef class P1Sampler2D(ElementSampler):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -51,7 +51,7 @@
cdef class P1Sampler3D(ElementSampler):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -67,7 +67,7 @@
# This typedef defines a function pointer that defines the system
# of equations that will be solved by the NonlinearSolveSamplers.
-#
+#
# inputs:
# x - pointer to the mapped coordinate
# vertices - pointer to the element vertices
@@ -78,15 +78,15 @@
# fx - the result of solving the system, should be close to 0
# once it is converged.
#
-ctypedef void (*func_type)(double* fx,
- double* x,
- double* vertices,
+ctypedef void (*func_type)(double* fx,
+ double* x,
+ double* vertices,
double* phys_x) nogil
# This typedef defines a function pointer that defines the Jacobian
-# matrix used by the NonlinearSolveSampler3D. Subclasses needed to
+# matrix used by the NonlinearSolveSampler3D. Subclasses needed to
# define a Jacobian function in this form.
-#
+#
# inputs:
# x - pointer to the mapped coordinate
# vertices - pointer to the element vertices
@@ -98,18 +98,18 @@
# scol - the second column of the jacobian
# tcol - the third column of the jaocobian
#
-ctypedef void (*jac_type3D)(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
+ctypedef void (*jac_type3D)(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
double* phys_x) nogil
# This typedef defines a function pointer that defines the Jacobian
-# matrix used by the NonlinearSolveSampler2D. Subclasses needed to
+# matrix used by the NonlinearSolveSampler2D. Subclasses needed to
# define a Jacobian function in this form.
-#
+#
# inputs:
# x - pointer to the mapped coordinate
# vertices - pointer to the element vertices
@@ -132,19 +132,19 @@
cdef int dim
cdef int max_iter
cdef np.float64_t tolerance
- cdef func_type func
+ cdef func_type func
cdef jac_type3D jac
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
-
+
cdef class Q1Sampler3D(NonlinearSolveSampler3D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -161,7 +161,7 @@
cdef class W1Sampler3D(NonlinearSolveSampler3D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -179,7 +179,7 @@
cdef class S2Sampler3D(NonlinearSolveSampler3D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -199,19 +199,19 @@
cdef int dim
cdef int max_iter
cdef np.float64_t tolerance
- cdef func_type func
+ cdef func_type func
cdef jac_type2D jac
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
-
+
cdef class Q1Sampler2D(NonlinearSolveSampler2D):
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil
@@ -221,3 +221,31 @@
double* vals) nogil
cdef int check_inside(self, double* mapped_coord) nogil
+
+cdef class T2Sampler2D(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 Tet2Sampler3D(NonlinearSolveSampler3D):
+
+ 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
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -1,6 +1,6 @@
"""
This file contains coordinate mappings between physical coordinates and those
-defined on unit elements, as well as doing the corresponding intracell
+defined on unit elements, as well as doing the corresponding intracell
interpolation on finite element data.
@@ -25,7 +25,11 @@
Q1Function2D, \
Q1Jacobian2D, \
W1Function3D, \
- W1Jacobian3D
+ W1Jacobian3D, \
+ T2Function2D, \
+ T2Jacobian2D, \
+ Tet2Function3D, \
+ Tet2Jacobian3D
cdef extern from "platform_dep.h":
double fmax(double x, double y) nogil
@@ -33,8 +37,8 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef double determinant_3x3(double* col0,
- double* col1,
+cdef double determinant_3x3(double* col0,
+ double* col1,
double* col2) nogil:
return col0[0]*col1[1]*col2[2] - col0[0]*col1[2]*col2[1] - \
col0[1]*col1[0]*col2[2] + col0[1]*col1[2]*col2[0] + \
@@ -49,7 +53,7 @@
cdef int i
err = fabs(f[0])
for i in range(1, dim):
- err = fmax(err, fabs(f[i]))
+ err = fmax(err, fabs(f[i]))
return err
@@ -58,7 +62,7 @@
This is a base class for sampling the value of a finite element solution
at an arbitrary point inside a mesh element. In general, this will be done
- by transforming the requested physical coordinate into a mapped coordinate
+ by transforming the requested physical coordinate into a mapped coordinate
system, sampling the solution in mapped coordinates, and returning the result.
This is not to be used directly; use one of the subclasses instead.
@@ -71,11 +75,11 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef void map_real_to_unit(self,
- double* mapped_x,
+ double* mapped_x,
double* vertices,
double* physical_x) nogil:
pass
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -115,7 +119,8 @@
'''
This implements sampling inside a linear, triangular mesh element.
- This mapping is linear and can be inverted easily.
+ This mapping is linear and can be inverted easily. Note that this
+ implementation uses triangular (or barycentric) coordinates.
'''
@@ -127,43 +132,43 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void map_real_to_unit(self, double* mapped_x,
+ cdef void map_real_to_unit(self, double* mapped_x,
double* vertices, double* physical_x) nogil:
-
+
cdef double[3] col0
cdef double[3] col1
cdef double[3] col2
-
+
col0[0] = vertices[0]
col0[1] = vertices[1]
col0[2] = 1.0
-
+
col1[0] = vertices[2]
col1[1] = vertices[3]
col1[2] = 1.0
-
+
col2[0] = vertices[4]
col2[1] = vertices[5]
col2[2] = 1.0
-
+
det = determinant_3x3(col0, col1, col2)
-
+
mapped_x[0] = ((vertices[3] - vertices[5])*physical_x[0] + \
(vertices[4] - vertices[2])*physical_x[1] + \
(vertices[2]*vertices[5] - vertices[4]*vertices[3])) / det
-
+
mapped_x[1] = ((vertices[5] - vertices[1])*physical_x[0] + \
(vertices[0] - vertices[4])*physical_x[1] + \
(vertices[4]*vertices[1] - vertices[0]*vertices[5])) / det
-
+
mapped_x[2] = 1.0 - mapped_x[1] - mapped_x[0]
-
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef double sample_at_unit_point(self,
- double* coord,
+ double* coord,
double* vals) nogil:
return vals[0]*coord[0] + vals[1]*coord[1] + vals[2]*coord[2]
@@ -197,16 +202,16 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- cdef void map_real_to_unit(self, double* mapped_x,
+ cdef void map_real_to_unit(self, double* mapped_x,
double* vertices, double* physical_x) nogil:
-
+
cdef int i
cdef double d
cdef double[3] bvec
cdef double[3] col0
cdef double[3] col1
cdef double[3] col2
-
+
# here, we express positions relative to the 4th element,
# which is selected by vertices[9]
for i in range(3):
@@ -214,7 +219,7 @@
col0[i] = vertices[0 + i] - vertices[9 + i]
col1[i] = vertices[3 + i] - vertices[9 + i]
col2[i] = vertices[6 + i] - vertices[9 + i]
-
+
d = determinant_3x3(col0, col1, col2)
mapped_x[0] = determinant_3x3(bvec, col1, col2)/d
mapped_x[1] = determinant_3x3(col0, bvec, col2)/d
@@ -225,7 +230,7 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef double sample_at_unit_point(self,
- double* coord,
+ double* coord,
double* vals) nogil:
return vals[0]*coord[0] + vals[1]*coord[1] + \
vals[2]*coord[2] + vals[3]*coord[3]
@@ -265,11 +270,11 @@
u = mapped_coord[1]
v = mapped_coord[2]
w = mapped_coord[0]
- if ((u < thresh) or
- (v < thresh) or
+ if ((u < thresh) or
+ (v < thresh) or
(w < thresh) or
- (fabs(u - 1) < thresh) or
- (fabs(v - 1) < thresh) or
+ (fabs(u - 1) < thresh) or
+ (fabs(v - 1) < thresh) or
(fabs(w - 1) < thresh)):
return 1
return -1
@@ -281,7 +286,7 @@
This is a base class for handling element samplers that require
a nonlinear solve to invert the mapping between coordinate systems.
- To do this, we perform Newton-Raphson iteration using a specified
+ To do this, we perform Newton-Raphson iteration using a specified
system of equations with an analytic Jacobian matrix. This solver
is hard-coded for 3D for reasons of efficiency. This is not to be
used directly, use one of the subclasses instead.
@@ -302,7 +307,7 @@
double* physical_x) nogil:
cdef int i
cdef double d
- cdef double[3] f
+ cdef double[3] f
cdef double[3] r
cdef double[3] s
cdef double[3] t
@@ -313,11 +318,11 @@
# initial guess
for i in range(3):
x[i] = 0.0
-
+
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 3)
-
+ err = maxnorm(f, 3)
+
# begin Newton iteration
while (err > self.tolerance and iterations < self.max_iter):
self.jac(r, s, t, x, vertices, physical_x)
@@ -325,7 +330,7 @@
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)
+ self.func(f, x, vertices, physical_x)
err = maxnorm(f, 3)
iterations += 1
@@ -340,7 +345,7 @@
cdef class Q1Sampler3D(NonlinearSolveSampler3D):
- '''
+ '''
This implements sampling inside a 3D, linear, hexahedral mesh element.
@@ -358,14 +363,14 @@
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
cdef double F, rm, rp, sm, sp, tm, tp
-
+
rm = 1.0 - coord[0]
rp = 1.0 + coord[0]
sm = 1.0 - coord[1]
sp = 1.0 + coord[1]
tm = 1.0 - coord[2]
tp = 1.0 + coord[2]
-
+
F = vals[0]*rm*sm*tm + vals[1]*rp*sm*tm + vals[2]*rp*sp*tm + vals[3]*rm*sp*tm + \
vals[4]*rm*sm*tp + vals[5]*rp*sm*tp + vals[6]*rp*sp*tp + vals[7]*rm*sp*tp
return 0.125*F
@@ -374,10 +379,10 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int check_inside(self, double* mapped_coord) nogil:
- # for hexes, the mapped coordinates all go from -1 to 1
+ # for hexes, the mapped coordinates all go from -1 to 1
# if we are inside the element.
if (fabs(mapped_coord[0]) - 1.0 > self.inclusion_tol or
- fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
+ fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
fabs(mapped_coord[2]) - 1.0 > self.inclusion_tol):
return 0
return 1
@@ -401,7 +406,7 @@
cdef class S2Sampler3D(NonlinearSolveSampler3D):
- '''
+ '''
This implements sampling inside a 3D, 20-node hexahedral mesh element.
@@ -459,7 +464,7 @@
@cython.cdivision(True)
cdef int check_inside(self, double* mapped_coord) nogil:
if (fabs(mapped_coord[0]) - 1.0 > self.inclusion_tol or
- fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
+ fabs(mapped_coord[1]) - 1.0 > self.inclusion_tol or
fabs(mapped_coord[2]) - 1.0 > self.inclusion_tol):
return 0
return 1
@@ -485,8 +490,8 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef inline void S2Function3D(double* fx,
- double* x,
- double* vertices,
+ double* x,
+ double* vertices,
double* phys_x) nogil:
cdef int i
cdef double r, s, t, rm, rp, sm, sp, tm, tp
@@ -533,9 +538,9 @@
cdef inline void S2Jacobian3D(double* rcol,
double* scol,
double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
cdef int i
cdef double r, s, t, rm, rp, sm, sp, tm, tp
@@ -616,7 +621,7 @@
cdef class W1Sampler3D(NonlinearSolveSampler3D):
- '''
+ '''
This implements sampling inside a 3D, linear, wedge mesh element.
@@ -648,15 +653,15 @@
@cython.wraparound(False)
@cython.cdivision(True)
cdef int check_inside(self, double* mapped_coord) nogil:
- # for wedges the bounds of the mapped coordinates are:
+ # for wedges the bounds of the mapped coordinates are:
# 0 <= mapped_coord[0] <= 1 - mapped_coord[1]
# 0 <= mapped_coord[1]
# -1 <= mapped_coord[2] <= 1
if (mapped_coord[0] < -self.inclusion_tol or
mapped_coord[0] + mapped_coord[1] - 1.0 > self.inclusion_tol):
- return 0
+ return 0
if (mapped_coord[1] < -self.inclusion_tol):
- return 0
+ return 0
if (fabs(mapped_coord[2]) - 1.0 > self.inclusion_tol):
return 0
return 1
@@ -675,7 +680,7 @@
near_edge_r = (r < thresh) or (fabs(r + s - 1.0) < thresh)
near_edge_s = (s < thresh)
near_edge_t = fabs(fabs(mapped_coord[2]) - 1.0) < thresh
-
+
# we use ray.instID to pass back whether the ray is near the
# element boundary or not (used to annotate mesh lines)
if (near_edge_r and near_edge_s):
@@ -694,7 +699,7 @@
This is a base class for handling element samplers that require
a nonlinear solve to invert the mapping between coordinate systems.
- To do this, we perform Newton-Raphson iteration using a specified
+ To do this, we perform Newton-Raphson iteration using a specified
system of equations with an analytic Jacobian matrix. This solver
is hard-coded for 2D for reasons of efficiency. This is not to be
used directly, use one of the subclasses instead.
@@ -724,20 +729,20 @@
# initial guess
for i in range(2):
x[i] = 0.0
-
+
# initial error norm
self.func(f, x, vertices, physical_x)
- err = maxnorm(f, 2)
-
+ err = maxnorm(f, 2)
+
# begin Newton iteration
while (err > 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
- self.func(f, x, vertices, physical_x)
+ self.func(f, x, vertices, physical_x)
err = maxnorm(f, 2)
iterations += 1
@@ -752,7 +757,7 @@
cdef class Q1Sampler2D(NonlinearSolveSampler2D):
- '''
+ '''
This implements sampling inside a 2D, linear, quadrilateral mesh element.
@@ -770,12 +775,12 @@
@cython.cdivision(True)
cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
cdef double F, rm, rp, sm, sp
-
+
rm = 1.0 - coord[0]
rp = 1.0 + coord[0]
sm = 1.0 - coord[1]
sp = 1.0 + coord[1]
-
+
F = vals[0]*rm*sm + vals[1]*rp*sm + vals[2]*rp*sp + vals[3]*rm*sp
return 0.25*F
@@ -790,6 +795,110 @@
return 0
return 1
+cdef class T2Sampler2D(NonlinearSolveSampler2D):
+
+ '''
+
+ This implements sampling inside a 2D, quadratic, triangular mesh
+ element. Note that this implementation uses canonical coordinates.
+
+ '''
+
+ def __init__(self):
+ super(T2Sampler2D, self).__init__()
+ self.num_mapped_coords = 2
+ self.dim = 2
+ self.func = T2Function2D
+ self.jac = T2Jacobian2D
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
+ cdef double phi0, phi1, phi2, phi3, phi4, phi5
+
+ phi0 = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
+ 2 * coord[1]**2 + 4 * coord[0] * coord[1]
+ phi1 = -coord[0] + 2 * coord[0]**2
+ phi2 = -coord[1] + 2 * coord[1]**2
+ phi3 = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1]
+ phi4 = 4 * coord[0] * coord[1]
+ phi5 = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1]
+
+ return vals[0]*phi0 + vals[1]*phi1 + vals[2]*phi2 + vals[3]*phi3 + \
+ vals[4]*phi4 + vals[5]*phi5
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_inside(self, double* mapped_coord) nogil:
+ # for canonical tris, we check whether the mapped_coords are between
+ # 0 and 1.
+ cdef int i
+ for i in range(2):
+ if (mapped_coord[i] < -self.inclusion_tol or
+ mapped_coord[i] - 1.0 > self.inclusion_tol):
+ return 0
+ return 1
+
+cdef class Tet2Sampler3D(NonlinearSolveSampler3D):
+
+ '''
+
+ This implements sampling inside a 3D, quadratic, tetrahedral mesh
+ element. Note that this implementation uses canonical coordinates.
+
+ '''
+
+ def __init__(self):
+ super(Tet2Sampler3D, self).__init__()
+ self.num_mapped_coords = 3
+ self.dim = 3
+ self.func = Tet2Function3D
+ self.jac = Tet2Jacobian3D
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef double sample_at_unit_point(self, double* coord, double* vals) nogil:
+ cdef double[10] phi
+ cdef int i
+ cdef double return_value = 0
+
+ phi[0] = 1 - 3 * coord[0] + 2 * coord[0]**2 - 3 * coord[1] + \
+ 2 * coord[1]**2 - 3 * coord[2] + 2 * coord[2]**2 + \
+ 4 * coord[0] * coord[1] + 4 * coord[0] * coord[2] + \
+ 4 * coord[1] * coord[2]
+ phi[1] = -coord[0] + 2 * coord[0]**2
+ phi[2] = -coord[1] + 2 * coord[1]**2
+ phi[3] = -coord[2] + 2 * coord[2]**2
+ phi[4] = 4 * coord[0] - 4 * coord[0]**2 - 4 * coord[0] * coord[1] - \
+ 4 * coord[0] * coord[2]
+ phi[5] = 4 * coord[0] * coord[1]
+ phi[6] = 4 * coord[1] - 4 * coord[1]**2 - 4 * coord[0] * coord[1] - \
+ 4 * coord[1] * coord[2]
+ phi[7] = 4 * coord[2] - 4 * coord[2]**2 - 4 * coord[2] * coord[0] - \
+ 4 * coord[2] * coord[1]
+ phi[8] = 4 * coord[0] * coord[2]
+ phi[9] = 4 * coord[1] * coord[2]
+
+ for i in range(10):
+ return_value += phi[i] * vals[i]
+
+ return return_value
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ @cython.cdivision(True)
+ cdef int check_inside(self, double* mapped_coord) nogil:
+ # for canonical tets, we check whether the mapped_coords are between
+ # 0 and 1.
+ cdef int i
+ for i in range(3):
+ if (mapped_coord[i] < -self.inclusion_tol or
+ mapped_coord[i] - 1.0 > self.inclusion_tol):
+ return 0
+ return 1
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -894,3 +1003,37 @@
<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):
+
+ cdef double val
+
+ sampler = T2Sampler2D()
+
+ 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_tet2_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 = Tet2Sampler3D()
+
+ val = sampler.sample_at_real_point(<double*> vertices.data,
+ <double*> field_values.data,
+ <double*> physical_x.data)
+
+ return val
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -30,7 +30,9 @@
S2Sampler3D, \
P1Sampler2D, \
Q1Sampler2D, \
- W1Sampler3D
+ W1Sampler3D, \
+ T2Sampler2D, \
+ Tet2Sampler3D
cdef extern from "pixelization_constants.h":
enum:
@@ -201,7 +203,7 @@
my_array[j,i] += (dsp * overlap1) * overlap2
else:
my_array[j,i] = dsp
-
+
return my_array
@cython.cdivision(True)
@@ -309,10 +311,10 @@
cdef np.float64_t r_i, theta_i, dr_i, dtheta_i, dthetamin
cdef np.float64_t costheta, sintheta
cdef int i, pi, pj
-
+
imax = radius.argmax()
rmax = radius[imax] + dradius[imax]
-
+
if input_img is None:
img = np.zeros((buff_size[0], buff_size[1]))
img[:] = np.nan
@@ -363,7 +365,7 @@
sintheta = math.sin(theta_i)
while r_i < r0 + dr_i:
if rmax <= r_i:
- r_i += 0.5*dx
+ r_i += 0.5*dx
continue
y = r_i * costheta
x = r_i * sintheta
@@ -374,7 +376,7 @@
if img[pi, pj] != img[pi, pj]:
img[pi, pj] = 0.0
img[pi, pj] = field[i]
- r_i += 0.5*dx
+ r_i += 0.5*dx
theta_i += dthetamin
return img
@@ -412,7 +414,7 @@
cdef np.float64_t s2 = math.sqrt(2.0)
cdef np.float64_t xmax, ymax, xmin, ymin
nf = field.shape[0]
-
+
if input_img is None:
img = np.zeros((buff_size[0], buff_size[1]))
img[:] = np.nan
@@ -551,7 +553,7 @@
np.ndarray[np.int64_t, ndim=2] conn,
buff_size,
np.ndarray[np.float64_t, ndim=2] field,
- extents,
+ extents,
int index_offset = 0):
cdef np.ndarray[np.float64_t, ndim=3] img
img = np.zeros(buff_size, dtype="float64")
@@ -594,6 +596,10 @@
sampler = P1Sampler2D()
elif ndim == 2 and nvertices == 4:
sampler = Q1Sampler2D()
+ elif ndim == 2 and nvertices == 6:
+ sampler = T2Sampler2D()
+ elif ndim == 3 and nvertices == 10:
+ sampler = Tet2Sampler3D()
else:
raise YTElementTypeNotRecognized(ndim, nvertices)
@@ -602,7 +608,7 @@
if buff_size[2] != 1:
raise RuntimeError("Slices of 2D datasets must be "
"perpendicular to the 'z' direction.")
-
+
# allocate temporary storage
vertices = <np.float64_t *> malloc(ndim * sizeof(np.float64_t) * nvertices)
field_vals = <np.float64_t *> malloc(sizeof(np.float64_t) * num_field_vals)
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac 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
@@ -22,7 +22,9 @@
test_tri_sampler, \
test_quad_sampler, \
test_hex20_sampler, \
- test_wedge_sampler
+ test_wedge_sampler, \
+ test_tri2_sampler, \
+ test_tet2_sampler
def check_all_vertices(sampler, vertices, field_values):
@@ -127,3 +129,34 @@
field_values = np.array([1., 2., 3., 4., 5., 6.])
check_all_vertices(test_wedge_sampler, vertices, field_values)
+
+def test_T2Sampler2D():
+
+ vertices = np.array([[0.1 , 0.2 ],
+ [0.3 , 0.5 ],
+ [0.2 , 0.9 ],
+ [0.2 , 0.35],
+ [0.25, 0.7 ],
+ [0.15, 0.55]])
+
+ field_values = np.array([15., 37., 49., 32., 46., 24.])
+
+ check_all_vertices(test_tri2_sampler, vertices, field_values)
+
+
+def test_Tet2Sampler3D():
+
+ vertices = np.array([[0.3 , -0.4 , 0.6] ,
+ [1.7 , -0.7 , 0.8] ,
+ [0.4 , 1.2 , 0.4] ,
+ [0.4 , -0.2 , 2.0] ,
+ [1.0 , -0.55 , 0.7] ,
+ [1.05 , 0.25 , 0.6] ,
+ [0.35 , 0.4 , 0.5] ,
+ [0.35 , -0.3 , 1.3] ,
+ [1.05 , -0.45 , 1.4] ,
+ [0.4 , 0.5 , 1.2]])
+
+ field_values = np.array([15., 37., 49., 24., 30., 44., 20., 17., 32., 36.])
+
+ check_all_vertices(test_tet2_sampler, vertices, field_values)
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -5,9 +5,9 @@
Usage (from the yt/utilities directory):
-python mesh_code_generation.py
+python mesh_code_generation.py
-This will generate the necessary functions and write them to
+This will generate the necessary functions and write them to
yt/utilities/lib/autogenerated_element_samplers.pyx.
@@ -120,7 +120,7 @@
self.jacobian_name = '%sJacobian%dD' % (self.mesh_type, self.num_dim)
if (self.num_dim == 3):
- self.jacobian_header = jac_def_template_3D % self.jacobian_name
+ self.jacobian_header = jac_def_template_3D % self.jacobian_name
self.jacobian_declaration = jac_dec_template_3D % self.jacobian_name
elif (self.num_dim == 2):
@@ -137,15 +137,15 @@
function_code = self.function_header
for i in range(self.num_dim):
function_code += '\t' + ccode(self.f[i, 0], self.fx[i, 0]) + '\n'
-
+
jacobian_code = self.jacobian_header
for i in range(self.num_dim):
jacobian_code += '\t' + ccode(self.J[i,0], self.rcol[i, 0]) + '\n'
jacobian_code += '\t' + ccode(self.J[i,1], self.scol[i, 0]) + '\n'
- if self.num_dim == 2:
+ if self.num_dim == 2:
continue
jacobian_code += '\t' + ccode(self.J[i,2], self.tcol[i, 0]) + '\n'
-
+
return function_code, jacobian_code
def get_interpolator_declaration(self):
@@ -169,9 +169,10 @@
pyx_file.write(file_header)
pyx_file.write("\n \n")
- pyx_file.write("cimport cython \n \n")
+ pyx_file.write("cimport cython \n")
+ pyx_file.write("from libc.math cimport pow \n")
pyx_file.write("\n \n")
-
+
for _, mesh_data in sorted(mesh_types.items()):
codegen = MeshCodeGenerator(mesh_data)
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -12,7 +12,7 @@
(1 + x[0])*(1 - x[1])*(1 + x[2])/8.,
(1 + x[0])*(1 + x[1])*(1 + x[2])/8.,
(1 - x[0])*(1 + x[1])*(1 + x[2])/8.]
-
+
Quad4:
mesh_type: Q1
num_dim: 2
@@ -35,4 +35,34 @@
x[1] * (1 - x[2]) / 2.,
(1 - x[0] - x[1]) * (1 + x[2]) / 2.,
x[0] * (1 + x[2]) / 2.,
- x[1] * (1 + x[2]) / 2.]
\ No newline at end of file
+ x[1] * (1 + x[2]) / 2.]
+
+Tri6:
+ mesh_type: T2
+ num_dim: 2
+ num_vertices: 6
+ num_mapped_coords: 2
+ shape_functions: |
+ [1 - 3 * x[0] + 2 * x[0]**2 - 3 * x[1] + 2 * x[1]**2 + 4 * x[0] * x[1],
+ -x[0] + 2 * x[0]**2,
+ -x[1] + 2 * x[1]**2,
+ 4 * x[0] - 4 * x[0]**2 - 4 * x[0] * x[1],
+ 4 * x[0] * x[1],
+ 4 * x[1] - 4 * x[1]**2 - 4 * x[0] * x[1]]
+
+Tet10:
+ mesh_type: Tet2
+ num_dim: 3
+ num_vertices: 10
+ num_mapped_coords: 3
+ shape_functions: |
+ [1 - 3 * x[0] + 2 * x[0]**2 - 3 * x[1] + 2 * x[1]**2 - 3 * x[2] + 2 * x[2]**2 + 4 * x[0] * x[1] + 4 * x[0] * x[2] + 4 * x[1] * x[2],
+ -x[0] + 2 * x[0]**2,
+ -x[1] + 2 * x[1]**2,
+ -x[2] + 2 * x[2]**2,
+ 4 * x[0] - 4 * x[0]**2 - 4 * x[0] * x[1] - 4 * x[0] * x[2],
+ 4 * x[0] * x[1],
+ 4 * x[1] - 4 * x[1]**2 - 4 * x[1] * x[0] - 4 * x[1] * x[2],
+ 4 * x[2] - 4 * x[2]**2 - 4 * x[2] * x[0] - 4 * x[2] * x[1],
+ 4 * x[0] * x[2],
+ 4 * x[1] * x[2]]
diff -r 3c53eabd6afdccb7bfc120884ff5014dfa599a2d -r b1e7a7e87d3fb46fb9e6ff90d762a7b258b6c9ac yt/visualization/tests/test_mesh_slices.py
--- a/yt/visualization/tests/test_mesh_slices.py
+++ b/yt/visualization/tests/test_mesh_slices.py
@@ -18,6 +18,10 @@
import yt
from yt.testing import fake_tetrahedral_ds
from yt.testing import fake_hexahedral_ds
+from yt.utilities.answer_testing.framework import \
+ requires_ds, \
+ data_dir_load, \
+ GenericImageTest
def setup():
@@ -25,6 +29,25 @@
from yt.config import ytcfg
ytcfg["yt", "__withintesting"] = "True"
+def compare(ds, field, test_prefix, decimals=12):
+ def slice_image(filename_prefix):
+ sl = yt.SlicePlot(ds, 'z', field)
+ sl.set_log('all', False)
+ image_file = sl.save(filename_prefix)
+ return image_file
+
+ slice_image.__name__ = "slice_{}".format(test_prefix)
+ test = GenericImageTest(ds, slice_image, decimals)
+ test.prefix = test_prefix
+ return test
+
+tri2 = "SecondOrderTris/RZ_p_no_parts_do_nothing_bcs_cone_out.e"
+
+ at requires_ds(tri2)
+def test_tri2():
+ ds = data_dir_load(tri2, kwargs={'step':-1})
+ for field in ds.field_list:
+ yield compare(ds, field, "answers_tri2_%s_%s" % (field[0], field[1]))
def test_mesh_slices():
# Perform I/O in safe place instead of yt main dir
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