[yt-svn] commit/yt: 19 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Aug 24 11:10:54 PDT 2016
19 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/831ebae13168/
Changeset: 831ebae13168
Branch: yt
User: atmyers
Date: 2016-08-02 20:58:23+00:00
Summary: adding the autogenerated sampler extension to setup.py
Affected #: 1 file
diff -r 2c73913367acd93482a88841daebdd6cece35d2c -r 831ebae13168af001422f00cd508ae07528f69ff setup.py
--- a/setup.py
+++ b/setup.py
@@ -114,6 +114,9 @@
["yt/utilities/spatial/ckdtree.pyx"],
include_dirs=["yt/utilities/lib/"],
libraries=std_libs),
+ Extension("yt.utilities.lib.autogenerated_element_samplers",
+ ["yt/utilities/lib/autogenerated_element_samplers.pyx"],
+ include_dirs=["yt/utilities/lib/"]),
Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs),
https://bitbucket.org/yt_analysis/yt/commits/04d403a17307/
Changeset: 04d403a17307
Branch: yt
User: atmyers
Date: 2016-08-02 20:58:53+00:00
Summary: Adding the class that does the codegen.
Affected #: 1 file
diff -r 831ebae13168af001422f00cd508ae07528f69ff -r 04d403a17307d3d329631f292892a2ff71f4558c yt/utilities/mesh_code_generation.py
--- /dev/null
+++ b/yt/utilities/mesh_code_generation.py
@@ -0,0 +1,150 @@
+from sympy import \
+ symarray, \
+ symbols, \
+ diff
+import re
+import yaml
+
+
+fun_dec_template = '''cdef inline void %s(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil \n'''
+
+fun_def_template = '''@cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void %s(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil: \n'''
+
+jac_dec_template = '''cdef inline void %s(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil \n'''
+
+jac_def_template = '''@cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void %s(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil: \n'''
+
+file_header = "# This file contains auto-generated functions for sampling \n" + \
+ "# inside finite element solutions for various mesh types."
+
+
+def replace_func(match):
+ s = match.group(0)
+ i = int(s[-3])
+ j = int(s[-1])
+ n = 3*j + i
+ return 'vertices[%d]' % n
+
+
+class MeshCodeGenerator:
+ def __init__(self, mesh_data):
+ self.mesh_type = mesh_data['mesh_type']
+ self.num_dim = mesh_data['num_dim']
+ self.num_vertices = mesh_data['num_vertices']
+ self.num_mapped_coords = mesh_data['num_mapped_coords']
+
+ x = symarray('x', self.num_mapped_coords)
+ self.x = x
+ self.N = eval(mesh_data['shape_functions'])
+ self._compute_jacobian()
+
+ def _compute_jacobian(self):
+
+ assert(self.num_vertices == len(self.N))
+ assert(self.num_dim == self.num_mapped_coords)
+
+ X = symarray('vertices', (self.num_dim, self.num_vertices))
+ physical_position = symbols(['phys_x[%s] ' % d for d in '012'[:self.num_dim]])
+
+ self.f = X.dot(self.N) - physical_position
+
+ self.J = symarray('J', (self.num_dim, self.num_dim))
+ for i in range(self.num_dim):
+ for j, var in enumerate(self.x):
+ self.J[i][j] = diff(self.f[i], var)
+
+ self.function_name = '%sFunction%dD' % (self.mesh_type, self.num_dim)
+ self.jacobian_name = '%sJacobian%dD' % (self.mesh_type, self.num_dim)
+
+ self.function_header = fun_def_template % self.function_name
+ self.jacobian_header = jac_def_template % self.jacobian_name
+
+ self.function_declaration = fun_dec_template % self.function_name
+ self.jacobian_declaration = jac_dec_template % self.jacobian_name
+
+ def get_function_line(self, i):
+ line = str(self.f[i])
+ for j in range(self.num_dim):
+ line = re.sub(r'x_%d' % j, 'x[%d]' % j, line)
+ line = re.sub(r'(vertices_._.)', replace_func, line)
+ return ''' fx[%d] = %s \n''' % (i, line)
+
+ def get_jacobian_line(self, i, j):
+ line = str(self.J[i, j])
+ for k in range(self.num_dim):
+ line = re.sub(r'x_%d' % k, 'x[%d]' % k, line)
+ line = re.sub(r'(vertices_._.)', replace_func, line)
+ col = 'rst'[j]
+ return ''' %scol[%d] = %s \n''' % (col, i, line)
+
+ def get_interpolator_definition(self):
+ function_code = self.function_header
+ for i in range(self.num_mapped_coords):
+ function_code += self.get_function_line(i)
+
+ jacobian_code = self.jacobian_header
+ for i in range(self.num_dim):
+ for j in range(self.num_dim):
+ jacobian_code += self.get_jacobian_line(i, j)
+
+ return function_code, jacobian_code
+
+ def get_interpolator_declaration(self):
+ return self.function_declaration, self.jacobian_declaration
+
+
+if __name__ == "__main__":
+
+ with open('mesh_types.yaml', 'r') as f:
+ lines = f.read()
+
+ mesh_types = yaml.load(lines)
+
+ pxd_file = open("lib/autogenerated_element_samplers.pxd", "w")
+ pyx_file = open("lib/autogenerated_element_samplers.pyx", "w")
+
+ pyx_file.write(file_header)
+ pyx_file.write("\n \n")
+ pyx_file.write("cimport cython \n \n")
+ pyx_file.write("\n \n")
+
+ for mesh_data in mesh_types.values():
+ codegen = MeshCodeGenerator(mesh_data)
+
+ function_code, jacobian_code = codegen.get_interpolator_definition()
+ function_decl, jacobian_decl = codegen.get_interpolator_declaration()
+
+ pxd_file.write(function_decl)
+ pxd_file.write("\n \n")
+ pxd_file.write(jacobian_decl)
+ pxd_file.write("\n \n")
+
+ pyx_file.write(function_code)
+ pyx_file.write("\n \n")
+ pyx_file.write(jacobian_code)
+ pyx_file.write("\n \n")
+
+ pxd_file.close()
+ pyx_file.close()
https://bitbucket.org/yt_analysis/yt/commits/b1d31f041750/
Changeset: b1d31f041750
Branch: yt
User: atmyers
Date: 2016-08-02 20:59:19+00:00
Summary: Adding the yaml file that contains the mesh metadata.
Affected #: 1 file
diff -r 04d403a17307d3d329631f292892a2ff71f4558c -r b1d31f0417508215a868f32f185081274350459e yt/utilities/mesh_types.yaml
--- /dev/null
+++ b/yt/utilities/mesh_types.yaml
@@ -0,0 +1,38 @@
+Hex8:
+ mesh_type: Q1
+ num_dim: 3
+ num_vertices: 8
+ num_mapped_coords: 3
+ shape_functions: |
+ [(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,
+ (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,
+ (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
+ num_vertices: 4
+ num_mapped_coords: 2
+ shape_functions: |
+ [(1 - x[0])*(1 - x[1])/4,
+ (1 + x[0])*(1 - x[1])/4,
+ (1 + x[0])*(1 + x[1])/4,
+ (1 - x[0])*(1 + x[1])/4]
+
+Wedge6:
+ mesh_type: W1
+ num_dim: 3
+ num_vertices: 6
+ num_mapped_coords: 3
+ shape_functions: |
+ [(1 - x[0] - x[2]) * (1 - x[2]) / 2,
+ x[0] * (1 - x[2]) / 2,
+ x[1] * (1 - x[2]) / 2,
+ (1 - x[0] - x[2]) * (1 + x[2]) / 2,
+ x[0] * (1 + x[2]) / 2,
+ x[1] * (1 + x[2]) / 2]
\ No newline at end of file
https://bitbucket.org/yt_analysis/yt/commits/3dd36f4187c4/
Changeset: 3dd36f4187c4
Branch: yt
User: atmyers
Date: 2016-08-02 20:59:42+00:00
Summary: checking in the autogenerated code itself.
Affected #: 2 files
diff -r b1d31f0417508215a868f32f185081274350459e -r 3dd36f4187c44d4406cfebd97228dda7b5fd30a3 yt/utilities/lib/autogenerated_element_samplers.pxd
--- /dev/null
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -0,0 +1,42 @@
+cdef inline void Q1Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef inline void Q1Jacobian2D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef inline void W1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef inline void W1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef inline void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef inline void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
diff -r b1d31f0417508215a868f32f185081274350459e -r 3dd36f4187c44d4406cfebd97228dda7b5fd30a3 yt/utilities/lib/autogenerated_element_samplers.pyx
--- /dev/null
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -0,0 +1,97 @@
+# This file contains auto-generated functions for sampling
+# inside finite element solutions for various mesh types.
+
+cimport cython
+
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void Q1Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = -phys_x[0] + vertices[0]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[3]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[6]*(x[0] + 1)*(x[1] + 1)/4 + vertices[9]*(-x[0] + 1)*(x[1] + 1)/4
+ fx[1] = -phys_x[1] + vertices[1]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[4]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[7]*(x[0] + 1)*(x[1] + 1)/4 + vertices[10]*(-x[0] + 1)*(x[1] + 1)/4
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void Q1Jacobian2D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -vertices[0]*(-x[1] + 1)/4 + vertices[3]*(-x[1] + 1)/4 + vertices[6]*(x[1] + 1)/4 - vertices[9]*(x[1] + 1)/4
+ scol[0] = -vertices[0]*(-x[0] + 1)/4 - vertices[3]*(x[0] + 1)/4 + vertices[6]*(x[0] + 1)/4 + vertices[9]*(-x[0] + 1)/4
+ rcol[1] = -vertices[1]*(-x[1] + 1)/4 + vertices[4]*(-x[1] + 1)/4 + vertices[7]*(x[1] + 1)/4 - vertices[10]*(x[1] + 1)/4
+ scol[1] = -vertices[1]*(-x[0] + 1)/4 - vertices[4]*(x[0] + 1)/4 + vertices[7]*(x[0] + 1)/4 + vertices[10]*(-x[0] + 1)/4
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void W1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = -phys_x[0] + vertices[0]*(-x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[3]*x[0]*(-x[2] + 1)/2 + vertices[6]*x[1]*(-x[2] + 1)/2 + vertices[9]*(x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[12]*x[0]*(x[2] + 1)/2 + vertices[15]*x[1]*(x[2] + 1)/2
+ fx[1] = -phys_x[1] + vertices[1]*(-x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[4]*x[0]*(-x[2] + 1)/2 + vertices[7]*x[1]*(-x[2] + 1)/2 + vertices[10]*(x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[13]*x[0]*(x[2] + 1)/2 + vertices[16]*x[1]*(x[2] + 1)/2
+ fx[2] = -phys_x[2] + vertices[2]*(-x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[5]*x[0]*(-x[2] + 1)/2 + vertices[8]*x[1]*(-x[2] + 1)/2 + vertices[11]*(x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[14]*x[0]*(x[2] + 1)/2 + vertices[17]*x[1]*(x[2] + 1)/2
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void W1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -vertices[0]*(-x[2] + 1)/2 + vertices[3]*(-x[2] + 1)/2 - vertices[9]*(x[2] + 1)/2 + vertices[12]*(x[2] + 1)/2
+ scol[0] = vertices[6]*(-x[2] + 1)/2 + vertices[15]*(x[2] + 1)/2
+ tcol[0] = -vertices[0]*(-x[2] + 1)/2 - vertices[0]*(-x[0] - x[2] + 1)/2 - vertices[3]*x[0]/2 - vertices[6]*x[1]/2 - vertices[9]*(x[2] + 1)/2 + vertices[9]*(-x[0] - x[2] + 1)/2 + vertices[12]*x[0]/2 + vertices[15]*x[1]/2
+ rcol[1] = -vertices[1]*(-x[2] + 1)/2 + vertices[4]*(-x[2] + 1)/2 - vertices[10]*(x[2] + 1)/2 + vertices[13]*(x[2] + 1)/2
+ scol[1] = vertices[7]*(-x[2] + 1)/2 + vertices[16]*(x[2] + 1)/2
+ tcol[1] = -vertices[1]*(-x[2] + 1)/2 - vertices[1]*(-x[0] - x[2] + 1)/2 - vertices[4]*x[0]/2 - vertices[7]*x[1]/2 - vertices[10]*(x[2] + 1)/2 + vertices[10]*(-x[0] - x[2] + 1)/2 + vertices[13]*x[0]/2 + vertices[16]*x[1]/2
+ rcol[2] = -vertices[2]*(-x[2] + 1)/2 + vertices[5]*(-x[2] + 1)/2 - vertices[11]*(x[2] + 1)/2 + vertices[14]*(x[2] + 1)/2
+ scol[2] = vertices[8]*(-x[2] + 1)/2 + vertices[17]*(x[2] + 1)/2
+ tcol[2] = -vertices[2]*(-x[2] + 1)/2 - vertices[2]*(-x[0] - x[2] + 1)/2 - vertices[5]*x[0]/2 - vertices[8]*x[1]/2 - vertices[11]*(x[2] + 1)/2 + vertices[11]*(-x[0] - x[2] + 1)/2 + vertices[14]*x[0]/2 + vertices[17]*x[1]/2
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = -phys_x[0] + vertices[0]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[3]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[9]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[12]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[15]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[18]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
+ fx[1] = -phys_x[1] + vertices[1]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[4]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[7]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[10]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[13]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[16]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[19]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[22]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
+ fx[2] = -phys_x[2] + vertices[2]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[5]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[8]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[11]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[14]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[17]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[20]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[23]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -vertices[0]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[3]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[9]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[12]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[15]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[18]*(x[1] + 1)*(x[2] + 1)/8 - vertices[21]*(x[1] + 1)*(x[2] + 1)/8
+ scol[0] = -vertices[0]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[3]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[9]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[12]*(-x[0] + 1)*(x[2] + 1)/8 - vertices[15]*(x[0] + 1)*(x[2] + 1)/8 + vertices[18]*(x[0] + 1)*(x[2] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[2] + 1)/8
+ tcol[0] = -vertices[0]*(-x[0] + 1)*(-x[1] + 1)/8 - vertices[3]*(x[0] + 1)*(-x[1] + 1)/8 - vertices[6]*(x[0] + 1)*(x[1] + 1)/8 - vertices[9]*(-x[0] + 1)*(x[1] + 1)/8 + vertices[12]*(-x[0] + 1)*(-x[1] + 1)/8 + vertices[15]*(x[0] + 1)*(-x[1] + 1)/8 + vertices[18]*(x[0] + 1)*(x[1] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[1] + 1)/8
+ rcol[1] = -vertices[1]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[4]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[7]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[10]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[13]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[16]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[19]*(x[1] + 1)*(x[2] + 1)/8 - vertices[22]*(x[1] + 1)*(x[2] + 1)/8
+ scol[1] = -vertices[1]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[4]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[7]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[10]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[13]*(-x[0] + 1)*(x[2] + 1)/8 - vertices[16]*(x[0] + 1)*(x[2] + 1)/8 + vertices[19]*(x[0] + 1)*(x[2] + 1)/8 + vertices[22]*(-x[0] + 1)*(x[2] + 1)/8
+ tcol[1] = -vertices[1]*(-x[0] + 1)*(-x[1] + 1)/8 - vertices[4]*(x[0] + 1)*(-x[1] + 1)/8 - vertices[7]*(x[0] + 1)*(x[1] + 1)/8 - vertices[10]*(-x[0] + 1)*(x[1] + 1)/8 + vertices[13]*(-x[0] + 1)*(-x[1] + 1)/8 + vertices[16]*(x[0] + 1)*(-x[1] + 1)/8 + vertices[19]*(x[0] + 1)*(x[1] + 1)/8 + vertices[22]*(-x[0] + 1)*(x[1] + 1)/8
+ rcol[2] = -vertices[2]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[5]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[8]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[11]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[14]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[17]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[20]*(x[1] + 1)*(x[2] + 1)/8 - vertices[23]*(x[1] + 1)*(x[2] + 1)/8
+ scol[2] = -vertices[2]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[5]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[8]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[11]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[14]*(-x[0] + 1)*(x[2] + 1)/8 - vertices[17]*(x[0] + 1)*(x[2] + 1)/8 + vertices[20]*(x[0] + 1)*(x[2] + 1)/8 + vertices[23]*(-x[0] + 1)*(x[2] + 1)/8
+ tcol[2] = -vertices[2]*(-x[0] + 1)*(-x[1] + 1)/8 - vertices[5]*(x[0] + 1)*(-x[1] + 1)/8 - vertices[8]*(x[0] + 1)*(x[1] + 1)/8 - vertices[11]*(-x[0] + 1)*(x[1] + 1)/8 + vertices[14]*(-x[0] + 1)*(-x[1] + 1)/8 + vertices[17]*(x[0] + 1)*(-x[1] + 1)/8 + vertices[20]*(x[0] + 1)*(x[1] + 1)/8 + vertices[23]*(-x[0] + 1)*(x[1] + 1)/8
+
+
https://bitbucket.org/yt_analysis/yt/commits/9553998e05c9/
Changeset: 9553998e05c9
Branch: yt
User: atmyers
Date: 2016-08-02 21:00:16+00:00
Summary: updating .hgignore
Affected #: 1 file
diff -r 3dd36f4187c44d4406cfebd97228dda7b5fd30a3 -r 9553998e05c994d02edc87599c27f255b025feb2 .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -23,6 +23,7 @@
yt/geometry/particle_smooth.c
yt/geometry/selection_routines.c
yt/utilities/amr_utils.c
+yt/utilities/lib/autogenerated_element_samplers.c
yt/utilities/kdtree/forthonf2c.h
yt/utilities/libconfig_wrapper.c
yt/utilities/spatial/ckdtree.c
https://bitbucket.org/yt_analysis/yt/commits/a84397180566/
Changeset: a84397180566
Branch: yt
User: atmyers
Date: 2016-08-02 21:44:07+00:00
Summary: some changes needed to get codegen to work in 2D.
Affected #: 1 file
diff -r 9553998e05c994d02edc87599c27f255b025feb2 -r a84397180566f80df66464f9c8c45d103115a6a3 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -19,14 +19,14 @@
double* vertices,
double* phys_x) nogil: \n'''
-jac_dec_template = '''cdef inline void %s(double* rcol,
+jac_dec_template_3D = '''cdef inline void %s(double* rcol,
double* scol,
double* tcol,
double* x,
double* vertices,
double* phys_x) nogil \n'''
-jac_def_template = '''@cython.boundscheck(False)
+jac_def_template_3D = '''@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef inline void %s(double* rcol,
@@ -36,18 +36,23 @@
double* vertices,
double* phys_x) nogil: \n'''
+jac_dec_template_2D = '''cdef inline void %s(double* A,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil \n'''
+
+jac_def_template_2D = '''@cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef inline void %s(double* A,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil: \n'''
+
file_header = "# This file contains auto-generated functions for sampling \n" + \
"# inside finite element solutions for various mesh types."
-def replace_func(match):
- s = match.group(0)
- i = int(s[-3])
- j = int(s[-1])
- n = 3*j + i
- return 'vertices[%d]' % n
-
-
class MeshCodeGenerator:
def __init__(self, mesh_data):
self.mesh_type = mesh_data['mesh_type']
@@ -76,28 +81,42 @@
self.J[i][j] = diff(self.f[i], var)
self.function_name = '%sFunction%dD' % (self.mesh_type, self.num_dim)
+ self.function_header = fun_def_template % self.function_name
+ self.function_declaration = fun_dec_template % self.function_name
+
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_declaration = jac_dec_template_3D % self.jacobian_name
+ elif (self.num_dim == 2):
+ self.jacobian_header = jac_def_template_2D % self.jacobian_name
+ self.jacobian_declaration = jac_dec_template_2D % self.jacobian_name
- self.function_header = fun_def_template % self.function_name
- self.jacobian_header = jac_def_template % self.jacobian_name
-
- self.function_declaration = fun_dec_template % self.function_name
- self.jacobian_declaration = jac_dec_template % self.jacobian_name
+ def replace_func(self, match):
+ s = match.group(0)
+ i = int(s[-3])
+ j = int(s[-1])
+ n = self.num_dim*j + i
+ return 'vertices[%d]' % n
def get_function_line(self, i):
line = str(self.f[i])
for j in range(self.num_dim):
line = re.sub(r'x_%d' % j, 'x[%d]' % j, line)
- line = re.sub(r'(vertices_._.)', replace_func, line)
+ line = re.sub(r'(vertices_._.)', self.replace_func, line)
return ''' fx[%d] = %s \n''' % (i, line)
def get_jacobian_line(self, i, j):
line = str(self.J[i, j])
for k in range(self.num_dim):
line = re.sub(r'x_%d' % k, 'x[%d]' % k, line)
- line = re.sub(r'(vertices_._.)', replace_func, line)
- col = 'rst'[j]
- return ''' %scol[%d] = %s \n''' % (col, i, line)
+ line = re.sub(r'(vertices_._.)', self.replace_func, line)
+ if (self.num_dim == 2):
+ return ''' A[%d] = %s \n''' % (2*i + j, line)
+ else:
+ assert(self.num_dim == 3)
+ col = 'rst'[j]
+ return ''' %scol[%d] = %s \n''' % (col, i, line)
def get_interpolator_definition(self):
function_code = self.function_header
https://bitbucket.org/yt_analysis/yt/commits/d936a0472a5a/
Changeset: d936a0472a5a
Branch: yt
User: atmyers
Date: 2016-08-02 22:01:01+00:00
Summary: Correcting the wedge shape functions.
Affected #: 1 file
diff -r a84397180566f80df66464f9c8c45d103115a6a3 -r d936a0472a5a867a3cdb256cc8931818319479de yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -30,9 +30,9 @@
num_vertices: 6
num_mapped_coords: 3
shape_functions: |
- [(1 - x[0] - x[2]) * (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,
- (1 - x[0] - x[2]) * (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
https://bitbucket.org/yt_analysis/yt/commits/b4fad5d8f7d0/
Changeset: b4fad5d8f7d0
Branch: yt
User: atmyers
Date: 2016-08-02 22:07:15+00:00
Summary: Removing the inline keyword to prevent compiler warnings in cythonized code.
Affected #: 1 file
diff -r d936a0472a5a867a3cdb256cc8931818319479de -r b4fad5d8f7d0bd912470ba08c6e1de7283c88416 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -6,48 +6,48 @@
import yaml
-fun_dec_template = '''cdef inline void %s(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil \n'''
+fun_dec_template = '''cdef void %s(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil \n'''
fun_def_template = '''@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void %s(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil: \n'''
+cdef void %s(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil: \n'''
-jac_dec_template_3D = '''cdef inline void %s(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil \n'''
+jac_dec_template_3D = '''cdef void %s(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil \n'''
jac_def_template_3D = '''@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void %s(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil: \n'''
+cdef void %s(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil: \n'''
-jac_dec_template_2D = '''cdef inline void %s(double* A,
- double* x,
- double* vertices,
- double* phys_x) nogil \n'''
+jac_dec_template_2D = '''cdef void %s(double* A,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil \n'''
jac_def_template_2D = '''@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void %s(double* A,
- double* x,
- double* vertices,
- double* phys_x) nogil: \n'''
+cdef void %s(double* A,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil: \n'''
file_header = "# This file contains auto-generated functions for sampling \n" + \
"# inside finite element solutions for various mesh types."
https://bitbucket.org/yt_analysis/yt/commits/8b013f7f8905/
Changeset: 8b013f7f8905
Branch: yt
User: atmyers
Date: 2016-08-02 22:08:04+00:00
Summary: Removing the old, hand-coded functions and jacobians from element_mappings.pyx
Affected #: 3 files
diff -r b4fad5d8f7d0bd912470ba08c6e1de7283c88416 -r 8b013f7f8905cc369fd97dbf4aa962b9b468811b yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -1,42 +1,40 @@
-cdef inline void Q1Function2D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil
+cdef void Q1Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
-cdef inline void Q1Jacobian2D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil
+cdef void Q1Jacobian2D(double* A,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
-cdef inline void W1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil
+cdef void W1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
-cdef inline void W1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- 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
-cdef inline void Q1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil
+cdef void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
-cdef inline void Q1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
diff -r b4fad5d8f7d0bd912470ba08c6e1de7283c88416 -r 8b013f7f8905cc369fd97dbf4aa962b9b468811b yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -8,68 +8,66 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void Q1Function2D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- fx[0] = -phys_x[0] + vertices[0]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[3]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[6]*(x[0] + 1)*(x[1] + 1)/4 + vertices[9]*(-x[0] + 1)*(x[1] + 1)/4
- fx[1] = -phys_x[1] + vertices[1]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[4]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[7]*(x[0] + 1)*(x[1] + 1)/4 + vertices[10]*(-x[0] + 1)*(x[1] + 1)/4
+cdef void Q1Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = -phys_x[0] + vertices[0]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[2]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[4]*(x[0] + 1)*(x[1] + 1)/4 + vertices[6]*(-x[0] + 1)*(x[1] + 1)/4
+ fx[1] = -phys_x[1] + vertices[1]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[3]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[5]*(x[0] + 1)*(x[1] + 1)/4 + vertices[7]*(-x[0] + 1)*(x[1] + 1)/4
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void Q1Jacobian2D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- rcol[0] = -vertices[0]*(-x[1] + 1)/4 + vertices[3]*(-x[1] + 1)/4 + vertices[6]*(x[1] + 1)/4 - vertices[9]*(x[1] + 1)/4
- scol[0] = -vertices[0]*(-x[0] + 1)/4 - vertices[3]*(x[0] + 1)/4 + vertices[6]*(x[0] + 1)/4 + vertices[9]*(-x[0] + 1)/4
- rcol[1] = -vertices[1]*(-x[1] + 1)/4 + vertices[4]*(-x[1] + 1)/4 + vertices[7]*(x[1] + 1)/4 - vertices[10]*(x[1] + 1)/4
- scol[1] = -vertices[1]*(-x[0] + 1)/4 - vertices[4]*(x[0] + 1)/4 + vertices[7]*(x[0] + 1)/4 + vertices[10]*(-x[0] + 1)/4
+cdef void Q1Jacobian2D(double* A,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ A[0] = -vertices[0]*(-x[1] + 1)/4 + vertices[2]*(-x[1] + 1)/4 + vertices[4]*(x[1] + 1)/4 - vertices[6]*(x[1] + 1)/4
+ A[1] = -vertices[0]*(-x[0] + 1)/4 - vertices[2]*(x[0] + 1)/4 + vertices[4]*(x[0] + 1)/4 + vertices[6]*(-x[0] + 1)/4
+ A[2] = -vertices[1]*(-x[1] + 1)/4 + vertices[3]*(-x[1] + 1)/4 + vertices[5]*(x[1] + 1)/4 - vertices[7]*(x[1] + 1)/4
+ A[3] = -vertices[1]*(-x[0] + 1)/4 - vertices[3]*(x[0] + 1)/4 + vertices[5]*(x[0] + 1)/4 + vertices[7]*(-x[0] + 1)/4
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void W1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- fx[0] = -phys_x[0] + vertices[0]*(-x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[3]*x[0]*(-x[2] + 1)/2 + vertices[6]*x[1]*(-x[2] + 1)/2 + vertices[9]*(x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[12]*x[0]*(x[2] + 1)/2 + vertices[15]*x[1]*(x[2] + 1)/2
- fx[1] = -phys_x[1] + vertices[1]*(-x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[4]*x[0]*(-x[2] + 1)/2 + vertices[7]*x[1]*(-x[2] + 1)/2 + vertices[10]*(x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[13]*x[0]*(x[2] + 1)/2 + vertices[16]*x[1]*(x[2] + 1)/2
- fx[2] = -phys_x[2] + vertices[2]*(-x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[5]*x[0]*(-x[2] + 1)/2 + vertices[8]*x[1]*(-x[2] + 1)/2 + vertices[11]*(x[2] + 1)*(-x[0] - x[2] + 1)/2 + vertices[14]*x[0]*(x[2] + 1)/2 + vertices[17]*x[1]*(x[2] + 1)/2
+cdef void W1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = -phys_x[0] + vertices[0]*(-x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[3]*x[0]*(-x[2] + 1)/2 + vertices[6]*x[1]*(-x[2] + 1)/2 + vertices[9]*(x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[12]*x[0]*(x[2] + 1)/2 + vertices[15]*x[1]*(x[2] + 1)/2
+ fx[1] = -phys_x[1] + vertices[1]*(-x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[4]*x[0]*(-x[2] + 1)/2 + vertices[7]*x[1]*(-x[2] + 1)/2 + vertices[10]*(x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[13]*x[0]*(x[2] + 1)/2 + vertices[16]*x[1]*(x[2] + 1)/2
+ fx[2] = -phys_x[2] + vertices[2]*(-x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[5]*x[0]*(-x[2] + 1)/2 + vertices[8]*x[1]*(-x[2] + 1)/2 + vertices[11]*(x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[14]*x[0]*(x[2] + 1)/2 + vertices[17]*x[1]*(x[2] + 1)/2
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void W1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- 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:
rcol[0] = -vertices[0]*(-x[2] + 1)/2 + vertices[3]*(-x[2] + 1)/2 - vertices[9]*(x[2] + 1)/2 + vertices[12]*(x[2] + 1)/2
- scol[0] = vertices[6]*(-x[2] + 1)/2 + vertices[15]*(x[2] + 1)/2
- tcol[0] = -vertices[0]*(-x[2] + 1)/2 - vertices[0]*(-x[0] - x[2] + 1)/2 - vertices[3]*x[0]/2 - vertices[6]*x[1]/2 - vertices[9]*(x[2] + 1)/2 + vertices[9]*(-x[0] - x[2] + 1)/2 + vertices[12]*x[0]/2 + vertices[15]*x[1]/2
+ scol[0] = -vertices[0]*(-x[2] + 1)/2 + vertices[6]*(-x[2] + 1)/2 - vertices[9]*(x[2] + 1)/2 + vertices[15]*(x[2] + 1)/2
+ tcol[0] = -vertices[0]*(-x[0] - x[1] + 1)/2 - vertices[3]*x[0]/2 - vertices[6]*x[1]/2 + vertices[9]*(-x[0] - x[1] + 1)/2 + vertices[12]*x[0]/2 + vertices[15]*x[1]/2
rcol[1] = -vertices[1]*(-x[2] + 1)/2 + vertices[4]*(-x[2] + 1)/2 - vertices[10]*(x[2] + 1)/2 + vertices[13]*(x[2] + 1)/2
- scol[1] = vertices[7]*(-x[2] + 1)/2 + vertices[16]*(x[2] + 1)/2
- tcol[1] = -vertices[1]*(-x[2] + 1)/2 - vertices[1]*(-x[0] - x[2] + 1)/2 - vertices[4]*x[0]/2 - vertices[7]*x[1]/2 - vertices[10]*(x[2] + 1)/2 + vertices[10]*(-x[0] - x[2] + 1)/2 + vertices[13]*x[0]/2 + vertices[16]*x[1]/2
+ scol[1] = -vertices[1]*(-x[2] + 1)/2 + vertices[7]*(-x[2] + 1)/2 - vertices[10]*(x[2] + 1)/2 + vertices[16]*(x[2] + 1)/2
+ tcol[1] = -vertices[1]*(-x[0] - x[1] + 1)/2 - vertices[4]*x[0]/2 - vertices[7]*x[1]/2 + vertices[10]*(-x[0] - x[1] + 1)/2 + vertices[13]*x[0]/2 + vertices[16]*x[1]/2
rcol[2] = -vertices[2]*(-x[2] + 1)/2 + vertices[5]*(-x[2] + 1)/2 - vertices[11]*(x[2] + 1)/2 + vertices[14]*(x[2] + 1)/2
- scol[2] = vertices[8]*(-x[2] + 1)/2 + vertices[17]*(x[2] + 1)/2
- tcol[2] = -vertices[2]*(-x[2] + 1)/2 - vertices[2]*(-x[0] - x[2] + 1)/2 - vertices[5]*x[0]/2 - vertices[8]*x[1]/2 - vertices[11]*(x[2] + 1)/2 + vertices[11]*(-x[0] - x[2] + 1)/2 + vertices[14]*x[0]/2 + vertices[17]*x[1]/2
+ scol[2] = -vertices[2]*(-x[2] + 1)/2 + vertices[8]*(-x[2] + 1)/2 - vertices[11]*(x[2] + 1)/2 + vertices[17]*(x[2] + 1)/2
+ tcol[2] = -vertices[2]*(-x[0] - x[1] + 1)/2 - vertices[5]*x[0]/2 - vertices[8]*x[1]/2 + vertices[11]*(-x[0] - x[1] + 1)/2 + vertices[14]*x[0]/2 + vertices[17]*x[1]/2
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void Q1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
+cdef void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
fx[0] = -phys_x[0] + vertices[0]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[3]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[9]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[12]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[15]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[18]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
fx[1] = -phys_x[1] + vertices[1]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[4]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[7]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[10]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[13]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[16]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[19]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[22]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
fx[2] = -phys_x[2] + vertices[2]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[5]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[8]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[11]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[14]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[17]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[20]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[23]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
@@ -78,12 +76,12 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef inline void Q1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
rcol[0] = -vertices[0]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[3]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[9]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[12]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[15]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[18]*(x[1] + 1)*(x[2] + 1)/8 - vertices[21]*(x[1] + 1)*(x[2] + 1)/8
scol[0] = -vertices[0]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[3]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[9]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[12]*(-x[0] + 1)*(x[2] + 1)/8 - vertices[15]*(x[0] + 1)*(x[2] + 1)/8 + vertices[18]*(x[0] + 1)*(x[2] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[2] + 1)/8
tcol[0] = -vertices[0]*(-x[0] + 1)*(-x[1] + 1)/8 - vertices[3]*(x[0] + 1)*(-x[1] + 1)/8 - vertices[6]*(x[0] + 1)*(x[1] + 1)/8 - vertices[9]*(-x[0] + 1)*(x[1] + 1)/8 + vertices[12]*(-x[0] + 1)*(-x[1] + 1)/8 + vertices[15]*(x[0] + 1)*(-x[1] + 1)/8 + vertices[18]*(x[0] + 1)*(x[1] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[1] + 1)/8
diff -r b4fad5d8f7d0bd912470ba08c6e1de7283c88416 -r 8b013f7f8905cc369fd97dbf4aa962b9b468811b yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -19,6 +19,13 @@
cimport cython
import numpy as np
from libc.math cimport fabs
+from yt.utilities.lib.autogenerated_element_samplers cimport \
+ Q1Function3D, \
+ Q1Jacobian3D, \
+ Q1Function2D, \
+ Q1Jacobian2D, \
+ W1Function3D, \
+ W1Jacobian3D
cdef extern from "platform_dep.h":
double fmax(double x, double y) nogil
@@ -392,70 +399,6 @@
return -1
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- cdef double rm, rp, sm, sp, tm, tp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
- tm = 1.0 - x[2]
- tp = 1.0 + x[2]
-
- for i in range(3):
- fx[i] = vertices[0 + i]*rm*sm*tm \
- + vertices[3 + i]*rp*sm*tm \
- + vertices[6 + i]*rp*sp*tm \
- + vertices[9 + i]*rm*sp*tm \
- + vertices[12 + i]*rm*sm*tp \
- + vertices[15 + i]*rp*sm*tp \
- + vertices[18 + i]*rp*sp*tp \
- + vertices[21 + i]*rm*sp*tp \
- - 8.0*phys_x[i]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- cdef double rm, rp, sm, sp, tm, tp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
- tm = 1.0 - x[2]
- tp = 1.0 + x[2]
-
- for i in range(3):
- rcol[i] = -sm*tm*vertices[0 + i] + sm*tm*vertices[3 + i] + \
- sp*tm*vertices[6 + i] - sp*tm*vertices[9 + i] - \
- sm*tp*vertices[12 + i] + sm*tp*vertices[15 + i] + \
- sp*tp*vertices[18 + i] - sp*tp*vertices[21 + i]
- scol[i] = -rm*tm*vertices[0 + i] - rp*tm*vertices[3 + i] + \
- rp*tm*vertices[6 + i] + rm*tm*vertices[9 + i] - \
- rm*tp*vertices[12 + i] - rp*tp*vertices[15 + i] + \
- rp*tp*vertices[18 + i] + rm*tp*vertices[21 + i]
- tcol[i] = -rm*sm*vertices[0 + i] - rp*sm*vertices[3 + i] - \
- rp*sp*vertices[6 + i] - rm*sp*vertices[9 + i] + \
- rm*sm*vertices[12 + i] + rp*sm*vertices[15 + i] + \
- rp*sp*vertices[18 + i] + rm*sp*vertices[21 + i]
-
-
cdef class S2Sampler3D(NonlinearSolveSampler3D):
'''
@@ -745,52 +688,6 @@
return -1
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void W1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- for i in range(3):
- fx[i] = vertices[0 + i]*(1.0 - x[0] - x[1])*(1.0 - x[2]) \
- + vertices[3 + i]*x[0]*(1.0 - x[2]) \
- + vertices[6 + i]*x[1]*(1.0 - x[2]) \
- + vertices[9 + i]*(1.0 - x[0] - x[1])*(1.0 + x[2]) \
- + vertices[12 + i]*x[0]*(1.0 + x[2]) \
- + vertices[15 + i]*x[1]*(1.0 + x[2]) \
- - 2.0*phys_x[i]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void W1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
-
- cdef int i
- for i in range(3):
- rcol[i] = (x[2] - 1.0) * vertices[0 + i] \
- - (x[2] - 1.0) * vertices[3 + i] \
- - (x[2] + 1.0) * vertices[9 + i] \
- + (x[2] + 1.0) * vertices[12 + i]
- scol[i] = (x[2] - 1.0) * vertices[0 + i] \
- - (x[2] - 1.0) * vertices[6 + i] \
- - (x[2] + 1.0) * vertices[9 + i] \
- + (x[2] + 1.0) * vertices[15 + i]
- tcol[i] = (x[0] + x[1] - 1.0) * vertices[0 + i] \
- - x[0] * vertices[3 + i] \
- - x[1] * vertices[6 + i] \
- - (x[0] + x[1] - 1.0) * vertices[9 + i] \
- + x[0] * vertices[12 + i] \
- + x[1] * vertices[15 + i]
-
-
cdef class NonlinearSolveSampler2D(ElementSampler):
'''
@@ -893,48 +790,6 @@
return 0
return 1
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Jacobian2D(double* A,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef double rm, rp, sm, sp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
-
- A[0] = -sm*vertices[0] + sm*vertices[2] + sp*vertices[4] - sp*vertices[6]
- A[1] = -rm*vertices[0] - rp*vertices[2] + rp*vertices[4] + rm*vertices[6]
- A[2] = -sm*vertices[1] + sm*vertices[3] + sp*vertices[5] - sp*vertices[7]
- A[3] = -rm*vertices[1] - rp*vertices[3] + rp*vertices[5] + rm*vertices[7]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Function2D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- cdef double rm, rp, sm, sp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
-
- for i in range(2):
- fx[i] = vertices[0 + i]*rm*sm \
- + vertices[2 + i]*rp*sm \
- + vertices[4 + i]*rp*sp \
- + vertices[6 + i]*rm*sp \
- - 4.0*phys_x[i]
-
@cython.boundscheck(False)
@cython.wraparound(False)
https://bitbucket.org/yt_analysis/yt/commits/e20f29c6f430/
Changeset: e20f29c6f430
Branch: yt
User: atmyers
Date: 2016-08-02 22:25:25+00:00
Summary: adding some docstrings.
Affected #: 2 files
diff -r 8b013f7f8905cc369fd97dbf4aa962b9b468811b -r e20f29c6f430d36f4f6277fbebf7059d369886b4 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -7,7 +7,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
+ at cython.cdivision(True)
cdef void Q1Function2D(double* fx,
double* x,
double* vertices,
@@ -18,7 +18,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
+ at cython.cdivision(True)
cdef void Q1Jacobian2D(double* A,
double* x,
double* vertices,
@@ -31,7 +31,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
+ at cython.cdivision(True)
cdef void W1Function3D(double* fx,
double* x,
double* vertices,
@@ -43,7 +43,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
+ at cython.cdivision(True)
cdef void W1Jacobian3D(double* rcol,
double* scol,
double* tcol,
@@ -63,7 +63,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
+ at cython.cdivision(True)
cdef void Q1Function3D(double* fx,
double* x,
double* vertices,
@@ -75,7 +75,7 @@
@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
+ at cython.cdivision(True)
cdef void Q1Jacobian3D(double* rcol,
double* scol,
double* tcol,
diff -r 8b013f7f8905cc369fd97dbf4aa962b9b468811b -r e20f29c6f430d36f4f6277fbebf7059d369886b4 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -1,3 +1,26 @@
+"""
+This file contains code for automatically generating the functions and jacobians
+used when sampling inside the supported finite element mesh types. The supported
+mesh types are defined in yt/utilities/mesh_types.yaml.
+
+Usage (from the yt/utilities directory:
+
+python mesh_code_generation.py
+
+This will generate the necessary functions and write them to
+yt/utilities/lib/autogenerated_element_samplers.pyx.
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
from sympy import \
symarray, \
symbols, \
@@ -6,55 +29,56 @@
import yaml
-fun_dec_template = '''cdef void %s(double* fx,
+# define some templates used below
+fun_signature = '''cdef void %s(double* fx,
double* x,
double* vertices,
- double* phys_x) nogil \n'''
+ double* phys_x) nogil'''
+fun_dec_template = fun_signature + ' \n'
fun_def_template = '''@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
-cdef void %s(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil: \n'''
+ at cython.cdivision(True) \n''' + fun_signature + ': \n'
-jac_dec_template_3D = '''cdef void %s(double* rcol,
+jac_signature_3D = '''cdef void %s(double* rcol,
double* scol,
double* tcol,
double* x,
double* vertices,
- double* phys_x) nogil \n'''
+ double* phys_x) nogil'''
+jac_dec_template_3D = jac_signature_3D + ' \n'
jac_def_template_3D = '''@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
-cdef void %s(double* rcol,
- double* scol,
- double* tcol,
+ at cython.cdivision(True) \n''' + jac_signature_3D + ': \n'
+
+jac_signature_2D = '''cdef void %s(double* A,
double* x,
double* vertices,
- double* phys_x) nogil: \n'''
-
-jac_dec_template_2D = '''cdef void %s(double* A,
- double* x,
- double* vertices,
- double* phys_x) nogil \n'''
-
+ double* phys_x) nogil'''
+jac_dec_template_2D = jac_signature_2D + ' \n'
jac_def_template_2D = '''@cython.boundscheck(False)
@cython.wraparound(False)
- at cython.cdivision(True)
-cdef void %s(double* A,
- double* x,
- double* vertices,
- double* phys_x) nogil: \n'''
+ at cython.cdivision(True) \n''' + jac_signature_2D + ': \n'
file_header = "# This file contains auto-generated functions for sampling \n" + \
"# inside finite element solutions for various mesh types."
class MeshCodeGenerator:
+ '''
+
+ A class for automatically generating the functions and jacobians used for
+ sampling inside finite element calculations.
+
+ '''
def __init__(self, mesh_data):
+ '''
+
+ Mesh data should be a dictionary containing information about the type
+ of elements used. See yt/utilities/mesh_types.yaml for more information.
+
+ '''
self.mesh_type = mesh_data['mesh_type']
self.num_dim = mesh_data['num_dim']
self.num_vertices = mesh_data['num_vertices']
@@ -92,25 +116,25 @@
self.jacobian_header = jac_def_template_2D % self.jacobian_name
self.jacobian_declaration = jac_dec_template_2D % self.jacobian_name
- def replace_func(self, match):
+ def _replace_func(self, match):
s = match.group(0)
i = int(s[-3])
j = int(s[-1])
n = self.num_dim*j + i
return 'vertices[%d]' % n
- def get_function_line(self, i):
+ def _get_function_line(self, i):
line = str(self.f[i])
for j in range(self.num_dim):
line = re.sub(r'x_%d' % j, 'x[%d]' % j, line)
- line = re.sub(r'(vertices_._.)', self.replace_func, line)
+ line = re.sub(r'(vertices_._.)', self._replace_func, line)
return ''' fx[%d] = %s \n''' % (i, line)
- def get_jacobian_line(self, i, j):
+ def _get_jacobian_line(self, i, j):
line = str(self.J[i, j])
for k in range(self.num_dim):
line = re.sub(r'x_%d' % k, 'x[%d]' % k, line)
- line = re.sub(r'(vertices_._.)', self.replace_func, line)
+ line = re.sub(r'(vertices_._.)', self._replace_func, line)
if (self.num_dim == 2):
return ''' A[%d] = %s \n''' % (2*i + j, line)
else:
@@ -119,18 +143,28 @@
return ''' %scol[%d] = %s \n''' % (col, i, line)
def get_interpolator_definition(self):
+ '''
+
+ This returns the function definitions for the given mesh type.
+
+ '''
function_code = self.function_header
for i in range(self.num_mapped_coords):
- function_code += self.get_function_line(i)
+ function_code += self._get_function_line(i)
jacobian_code = self.jacobian_header
for i in range(self.num_dim):
for j in range(self.num_dim):
- jacobian_code += self.get_jacobian_line(i, j)
+ jacobian_code += self._get_jacobian_line(i, j)
return function_code, jacobian_code
def get_interpolator_declaration(self):
+ '''
+
+ This returns the function declarations for the given mesh type.
+
+ '''
return self.function_declaration, self.jacobian_declaration
https://bitbucket.org/yt_analysis/yt/commits/5b7c5e5a3a22/
Changeset: 5b7c5e5a3a22
Branch: yt
User: atmyers
Date: 2016-08-02 22:34:23+00:00
Summary: Missed a ')'
Affected #: 1 file
diff -r e20f29c6f430d36f4f6277fbebf7059d369886b4 -r 5b7c5e5a3a22d4cd1dbbaca276ca007ac1c17458 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -3,7 +3,7 @@
used when sampling inside the supported finite element mesh types. The supported
mesh types are defined in yt/utilities/mesh_types.yaml.
-Usage (from the yt/utilities directory:
+Usage (from the yt/utilities directory):
python mesh_code_generation.py
https://bitbucket.org/yt_analysis/yt/commits/5b37f4313560/
Changeset: 5b37f4313560
Branch: yt
User: atmyers
Date: 2016-08-03 17:41:42+00:00
Summary: Adding mesh_types.yaml to MANIFEST.in
Affected #: 2 files
diff -r 5b7c5e5a3a22d4cd1dbbaca276ca007ac1c17458 -r 5b37f431356035e27bca10a03dd4a7b78df0c980 MANIFEST.in
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -3,6 +3,7 @@
include yt/visualization/mapserver/html/leaflet/*.css
include yt/visualization/mapserver/html/leaflet/*.js
include yt/visualization/mapserver/html/leaflet/images/*.png
+include yt/utilities/mesh_types.yaml
exclude scripts/pr_backport.py
recursive-include yt *.py *.pyx *.pxd *.h README* *.txt LICENSE* *.cu
recursive-include doc *.rst *.txt *.py *.ipynb *.png *.jpg *.css *.html
diff -r 5b7c5e5a3a22d4cd1dbbaca276ca007ac1c17458 -r 5b37f431356035e27bca10a03dd4a7b78df0c980 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -115,7 +115,6 @@
elif (self.num_dim == 2):
self.jacobian_header = jac_def_template_2D % self.jacobian_name
self.jacobian_declaration = jac_dec_template_2D % self.jacobian_name
-
def _replace_func(self, match):
s = match.group(0)
i = int(s[-3])
https://bitbucket.org/yt_analysis/yt/commits/3291782eb194/
Changeset: 3291782eb194
Branch: yt
User: atmyers
Date: 2016-08-03 17:57:23+00:00
Summary: Using ccode instead of str to print.
Affected #: 3 files
diff -r 5b37f431356035e27bca10a03dd4a7b78df0c980 -r 3291782eb19463320b8fdccba398fe66e98d87a1 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -12,8 +12,8 @@
double* x,
double* vertices,
double* phys_x) nogil:
- fx[0] = -phys_x[0] + vertices[0]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[2]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[4]*(x[0] + 1)*(x[1] + 1)/4 + vertices[6]*(-x[0] + 1)*(x[1] + 1)/4
- fx[1] = -phys_x[1] + vertices[1]*(-x[0] + 1)*(-x[1] + 1)/4 + vertices[3]*(x[0] + 1)*(-x[1] + 1)/4 + vertices[5]*(x[0] + 1)*(x[1] + 1)/4 + vertices[7]*(-x[0] + 1)*(x[1] + 1)/4
+ fx[0] = -phys_x[0] + 0.25*vertices[0]*(-x[0] + 1)*(-x[1] + 1) + 0.25*vertices[2]*(x[0] + 1)*(-x[1] + 1) + 0.25*vertices[4]*(x[0] + 1)*(x[1] + 1) + 0.25*vertices[6]*(-x[0] + 1)*(x[1] + 1)
+ fx[1] = -phys_x[1] + 0.25*vertices[1]*(-x[0] + 1)*(-x[1] + 1) + 0.25*vertices[3]*(x[0] + 1)*(-x[1] + 1) + 0.25*vertices[5]*(x[0] + 1)*(x[1] + 1) + 0.25*vertices[7]*(-x[0] + 1)*(x[1] + 1)
@cython.boundscheck(False)
@@ -23,10 +23,10 @@
double* x,
double* vertices,
double* phys_x) nogil:
- A[0] = -vertices[0]*(-x[1] + 1)/4 + vertices[2]*(-x[1] + 1)/4 + vertices[4]*(x[1] + 1)/4 - vertices[6]*(x[1] + 1)/4
- A[1] = -vertices[0]*(-x[0] + 1)/4 - vertices[2]*(x[0] + 1)/4 + vertices[4]*(x[0] + 1)/4 + vertices[6]*(-x[0] + 1)/4
- A[2] = -vertices[1]*(-x[1] + 1)/4 + vertices[3]*(-x[1] + 1)/4 + vertices[5]*(x[1] + 1)/4 - vertices[7]*(x[1] + 1)/4
- A[3] = -vertices[1]*(-x[0] + 1)/4 - vertices[3]*(x[0] + 1)/4 + vertices[5]*(x[0] + 1)/4 + vertices[7]*(-x[0] + 1)/4
+ A[0] = -0.25*vertices[0]*(-x[1] + 1) + 0.25*vertices[2]*(-x[1] + 1) + 0.25*vertices[4]*(x[1] + 1) - 0.25*vertices[6]*(x[1] + 1)
+ A[1] = -0.25*vertices[0]*(-x[0] + 1) - 0.25*vertices[2]*(x[0] + 1) + 0.25*vertices[4]*(x[0] + 1) + 0.25*vertices[6]*(-x[0] + 1)
+ A[2] = -0.25*vertices[1]*(-x[1] + 1) + 0.25*vertices[3]*(-x[1] + 1) + 0.25*vertices[5]*(x[1] + 1) - 0.25*vertices[7]*(x[1] + 1)
+ A[3] = -0.25*vertices[1]*(-x[0] + 1) - 0.25*vertices[3]*(x[0] + 1) + 0.25*vertices[5]*(x[0] + 1) + 0.25*vertices[7]*(-x[0] + 1)
@cython.boundscheck(False)
@@ -36,9 +36,9 @@
double* x,
double* vertices,
double* phys_x) nogil:
- fx[0] = -phys_x[0] + vertices[0]*(-x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[3]*x[0]*(-x[2] + 1)/2 + vertices[6]*x[1]*(-x[2] + 1)/2 + vertices[9]*(x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[12]*x[0]*(x[2] + 1)/2 + vertices[15]*x[1]*(x[2] + 1)/2
- fx[1] = -phys_x[1] + vertices[1]*(-x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[4]*x[0]*(-x[2] + 1)/2 + vertices[7]*x[1]*(-x[2] + 1)/2 + vertices[10]*(x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[13]*x[0]*(x[2] + 1)/2 + vertices[16]*x[1]*(x[2] + 1)/2
- fx[2] = -phys_x[2] + vertices[2]*(-x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[5]*x[0]*(-x[2] + 1)/2 + vertices[8]*x[1]*(-x[2] + 1)/2 + vertices[11]*(x[2] + 1)*(-x[0] - x[1] + 1)/2 + vertices[14]*x[0]*(x[2] + 1)/2 + vertices[17]*x[1]*(x[2] + 1)/2
+ fx[0] = -phys_x[0] + 0.5*vertices[0]*(-x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[3]*x[0]*(-x[2] + 1) + 0.5*vertices[6]*x[1]*(-x[2] + 1) + 0.5*vertices[9]*(x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[12]*x[0]*(x[2] + 1) + 0.5*vertices[15]*x[1]*(x[2] + 1)
+ fx[1] = -phys_x[1] + 0.5*vertices[1]*(-x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[4]*x[0]*(-x[2] + 1) + 0.5*vertices[7]*x[1]*(-x[2] + 1) + 0.5*vertices[10]*(x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[13]*x[0]*(x[2] + 1) + 0.5*vertices[16]*x[1]*(x[2] + 1)
+ fx[2] = -phys_x[2] + 0.5*vertices[2]*(-x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[5]*x[0]*(-x[2] + 1) + 0.5*vertices[8]*x[1]*(-x[2] + 1) + 0.5*vertices[11]*(x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[14]*x[0]*(x[2] + 1) + 0.5*vertices[17]*x[1]*(x[2] + 1)
@cython.boundscheck(False)
@@ -50,15 +50,15 @@
double* x,
double* vertices,
double* phys_x) nogil:
- rcol[0] = -vertices[0]*(-x[2] + 1)/2 + vertices[3]*(-x[2] + 1)/2 - vertices[9]*(x[2] + 1)/2 + vertices[12]*(x[2] + 1)/2
- scol[0] = -vertices[0]*(-x[2] + 1)/2 + vertices[6]*(-x[2] + 1)/2 - vertices[9]*(x[2] + 1)/2 + vertices[15]*(x[2] + 1)/2
- tcol[0] = -vertices[0]*(-x[0] - x[1] + 1)/2 - vertices[3]*x[0]/2 - vertices[6]*x[1]/2 + vertices[9]*(-x[0] - x[1] + 1)/2 + vertices[12]*x[0]/2 + vertices[15]*x[1]/2
- rcol[1] = -vertices[1]*(-x[2] + 1)/2 + vertices[4]*(-x[2] + 1)/2 - vertices[10]*(x[2] + 1)/2 + vertices[13]*(x[2] + 1)/2
- scol[1] = -vertices[1]*(-x[2] + 1)/2 + vertices[7]*(-x[2] + 1)/2 - vertices[10]*(x[2] + 1)/2 + vertices[16]*(x[2] + 1)/2
- tcol[1] = -vertices[1]*(-x[0] - x[1] + 1)/2 - vertices[4]*x[0]/2 - vertices[7]*x[1]/2 + vertices[10]*(-x[0] - x[1] + 1)/2 + vertices[13]*x[0]/2 + vertices[16]*x[1]/2
- rcol[2] = -vertices[2]*(-x[2] + 1)/2 + vertices[5]*(-x[2] + 1)/2 - vertices[11]*(x[2] + 1)/2 + vertices[14]*(x[2] + 1)/2
- scol[2] = -vertices[2]*(-x[2] + 1)/2 + vertices[8]*(-x[2] + 1)/2 - vertices[11]*(x[2] + 1)/2 + vertices[17]*(x[2] + 1)/2
- tcol[2] = -vertices[2]*(-x[0] - x[1] + 1)/2 - vertices[5]*x[0]/2 - vertices[8]*x[1]/2 + vertices[11]*(-x[0] - x[1] + 1)/2 + vertices[14]*x[0]/2 + vertices[17]*x[1]/2
+ rcol[0] = -0.5*vertices[0]*(-x[2] + 1) + 0.5*vertices[3]*(-x[2] + 1) - 0.5*vertices[9]*(x[2] + 1) + 0.5*vertices[12]*(x[2] + 1)
+ scol[0] = -0.5*vertices[0]*(-x[2] + 1) + 0.5*vertices[6]*(-x[2] + 1) - 0.5*vertices[9]*(x[2] + 1) + 0.5*vertices[15]*(x[2] + 1)
+ tcol[0] = -0.5*vertices[0]*(-x[0] - x[1] + 1) - 0.5*vertices[3]*x[0] - 0.5*vertices[6]*x[1] + 0.5*vertices[9]*(-x[0] - x[1] + 1) + 0.5*vertices[12]*x[0] + 0.5*vertices[15]*x[1]
+ rcol[1] = -0.5*vertices[1]*(-x[2] + 1) + 0.5*vertices[4]*(-x[2] + 1) - 0.5*vertices[10]*(x[2] + 1) + 0.5*vertices[13]*(x[2] + 1)
+ scol[1] = -0.5*vertices[1]*(-x[2] + 1) + 0.5*vertices[7]*(-x[2] + 1) - 0.5*vertices[10]*(x[2] + 1) + 0.5*vertices[16]*(x[2] + 1)
+ tcol[1] = -0.5*vertices[1]*(-x[0] - x[1] + 1) - 0.5*vertices[4]*x[0] - 0.5*vertices[7]*x[1] + 0.5*vertices[10]*(-x[0] - x[1] + 1) + 0.5*vertices[13]*x[0] + 0.5*vertices[16]*x[1]
+ rcol[2] = -0.5*vertices[2]*(-x[2] + 1) + 0.5*vertices[5]*(-x[2] + 1) - 0.5*vertices[11]*(x[2] + 1) + 0.5*vertices[14]*(x[2] + 1)
+ scol[2] = -0.5*vertices[2]*(-x[2] + 1) + 0.5*vertices[8]*(-x[2] + 1) - 0.5*vertices[11]*(x[2] + 1) + 0.5*vertices[17]*(x[2] + 1)
+ tcol[2] = -0.5*vertices[2]*(-x[0] - x[1] + 1) - 0.5*vertices[5]*x[0] - 0.5*vertices[8]*x[1] + 0.5*vertices[11]*(-x[0] - x[1] + 1) + 0.5*vertices[14]*x[0] + 0.5*vertices[17]*x[1]
@cython.boundscheck(False)
@@ -68,9 +68,9 @@
double* x,
double* vertices,
double* phys_x) nogil:
- fx[0] = -phys_x[0] + vertices[0]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[3]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[9]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[12]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[15]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[18]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
- fx[1] = -phys_x[1] + vertices[1]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[4]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[7]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[10]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[13]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[16]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[19]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[22]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
- fx[2] = -phys_x[2] + vertices[2]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[5]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[8]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[11]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1)/8 + vertices[14]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[17]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1)/8 + vertices[20]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8 + vertices[23]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)/8
+ fx[0] = -phys_x[0] + 0.125*vertices[0]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[3]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[6]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[9]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[12]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[15]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[18]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1) + 0.125*vertices[21]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)
+ fx[1] = -phys_x[1] + 0.125*vertices[1]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[4]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[7]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[10]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[13]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[16]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[19]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1) + 0.125*vertices[22]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)
+ fx[2] = -phys_x[2] + 0.125*vertices[2]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[5]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[8]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[11]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[14]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[17]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[20]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1) + 0.125*vertices[23]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)
@cython.boundscheck(False)
@@ -82,14 +82,14 @@
double* x,
double* vertices,
double* phys_x) nogil:
- rcol[0] = -vertices[0]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[3]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[9]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[12]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[15]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[18]*(x[1] + 1)*(x[2] + 1)/8 - vertices[21]*(x[1] + 1)*(x[2] + 1)/8
- scol[0] = -vertices[0]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[3]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[6]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[9]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[12]*(-x[0] + 1)*(x[2] + 1)/8 - vertices[15]*(x[0] + 1)*(x[2] + 1)/8 + vertices[18]*(x[0] + 1)*(x[2] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[2] + 1)/8
- tcol[0] = -vertices[0]*(-x[0] + 1)*(-x[1] + 1)/8 - vertices[3]*(x[0] + 1)*(-x[1] + 1)/8 - vertices[6]*(x[0] + 1)*(x[1] + 1)/8 - vertices[9]*(-x[0] + 1)*(x[1] + 1)/8 + vertices[12]*(-x[0] + 1)*(-x[1] + 1)/8 + vertices[15]*(x[0] + 1)*(-x[1] + 1)/8 + vertices[18]*(x[0] + 1)*(x[1] + 1)/8 + vertices[21]*(-x[0] + 1)*(x[1] + 1)/8
- rcol[1] = -vertices[1]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[4]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[7]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[10]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[13]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[16]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[19]*(x[1] + 1)*(x[2] + 1)/8 - vertices[22]*(x[1] + 1)*(x[2] + 1)/8
- scol[1] = -vertices[1]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[4]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[7]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[10]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[13]*(-x[0] + 1)*(x[2] + 1)/8 - vertices[16]*(x[0] + 1)*(x[2] + 1)/8 + vertices[19]*(x[0] + 1)*(x[2] + 1)/8 + vertices[22]*(-x[0] + 1)*(x[2] + 1)/8
- tcol[1] = -vertices[1]*(-x[0] + 1)*(-x[1] + 1)/8 - vertices[4]*(x[0] + 1)*(-x[1] + 1)/8 - vertices[7]*(x[0] + 1)*(x[1] + 1)/8 - vertices[10]*(-x[0] + 1)*(x[1] + 1)/8 + vertices[13]*(-x[0] + 1)*(-x[1] + 1)/8 + vertices[16]*(x[0] + 1)*(-x[1] + 1)/8 + vertices[19]*(x[0] + 1)*(x[1] + 1)/8 + vertices[22]*(-x[0] + 1)*(x[1] + 1)/8
- rcol[2] = -vertices[2]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[5]*(-x[1] + 1)*(-x[2] + 1)/8 + vertices[8]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[11]*(x[1] + 1)*(-x[2] + 1)/8 - vertices[14]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[17]*(-x[1] + 1)*(x[2] + 1)/8 + vertices[20]*(x[1] + 1)*(x[2] + 1)/8 - vertices[23]*(x[1] + 1)*(x[2] + 1)/8
- scol[2] = -vertices[2]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[5]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[8]*(x[0] + 1)*(-x[2] + 1)/8 + vertices[11]*(-x[0] + 1)*(-x[2] + 1)/8 - vertices[14]*(-x[0] + 1)*(x[2] + 1)/8 - vertices[17]*(x[0] + 1)*(x[2] + 1)/8 + vertices[20]*(x[0] + 1)*(x[2] + 1)/8 + vertices[23]*(-x[0] + 1)*(x[2] + 1)/8
- tcol[2] = -vertices[2]*(-x[0] + 1)*(-x[1] + 1)/8 - vertices[5]*(x[0] + 1)*(-x[1] + 1)/8 - vertices[8]*(x[0] + 1)*(x[1] + 1)/8 - vertices[11]*(-x[0] + 1)*(x[1] + 1)/8 + vertices[14]*(-x[0] + 1)*(-x[1] + 1)/8 + vertices[17]*(x[0] + 1)*(-x[1] + 1)/8 + vertices[20]*(x[0] + 1)*(x[1] + 1)/8 + vertices[23]*(-x[0] + 1)*(x[1] + 1)/8
+ rcol[0] = -0.125*vertices[0]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[3]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[6]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[9]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[12]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[15]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[18]*(x[1] + 1)*(x[2] + 1) - 0.125*vertices[21]*(x[1] + 1)*(x[2] + 1)
+ scol[0] = -0.125*vertices[0]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[3]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[6]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[9]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[12]*(-x[0] + 1)*(x[2] + 1) - 0.125*vertices[15]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[18]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[21]*(-x[0] + 1)*(x[2] + 1)
+ tcol[0] = -0.125*vertices[0]*(-x[0] + 1)*(-x[1] + 1) - 0.125*vertices[3]*(x[0] + 1)*(-x[1] + 1) - 0.125*vertices[6]*(x[0] + 1)*(x[1] + 1) - 0.125*vertices[9]*(-x[0] + 1)*(x[1] + 1) + 0.125*vertices[12]*(-x[0] + 1)*(-x[1] + 1) + 0.125*vertices[15]*(x[0] + 1)*(-x[1] + 1) + 0.125*vertices[18]*(x[0] + 1)*(x[1] + 1) + 0.125*vertices[21]*(-x[0] + 1)*(x[1] + 1)
+ rcol[1] = -0.125*vertices[1]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[4]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[7]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[10]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[13]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[16]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[19]*(x[1] + 1)*(x[2] + 1) - 0.125*vertices[22]*(x[1] + 1)*(x[2] + 1)
+ scol[1] = -0.125*vertices[1]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[4]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[7]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[10]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[13]*(-x[0] + 1)*(x[2] + 1) - 0.125*vertices[16]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[19]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[22]*(-x[0] + 1)*(x[2] + 1)
+ tcol[1] = -0.125*vertices[1]*(-x[0] + 1)*(-x[1] + 1) - 0.125*vertices[4]*(x[0] + 1)*(-x[1] + 1) - 0.125*vertices[7]*(x[0] + 1)*(x[1] + 1) - 0.125*vertices[10]*(-x[0] + 1)*(x[1] + 1) + 0.125*vertices[13]*(-x[0] + 1)*(-x[1] + 1) + 0.125*vertices[16]*(x[0] + 1)*(-x[1] + 1) + 0.125*vertices[19]*(x[0] + 1)*(x[1] + 1) + 0.125*vertices[22]*(-x[0] + 1)*(x[1] + 1)
+ rcol[2] = -0.125*vertices[2]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[5]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[8]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[11]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[14]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[17]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[20]*(x[1] + 1)*(x[2] + 1) - 0.125*vertices[23]*(x[1] + 1)*(x[2] + 1)
+ scol[2] = -0.125*vertices[2]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[5]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[8]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[11]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[14]*(-x[0] + 1)*(x[2] + 1) - 0.125*vertices[17]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[20]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[23]*(-x[0] + 1)*(x[2] + 1)
+ tcol[2] = -0.125*vertices[2]*(-x[0] + 1)*(-x[1] + 1) - 0.125*vertices[5]*(x[0] + 1)*(-x[1] + 1) - 0.125*vertices[8]*(x[0] + 1)*(x[1] + 1) - 0.125*vertices[11]*(-x[0] + 1)*(x[1] + 1) + 0.125*vertices[14]*(-x[0] + 1)*(-x[1] + 1) + 0.125*vertices[17]*(x[0] + 1)*(-x[1] + 1) + 0.125*vertices[20]*(x[0] + 1)*(x[1] + 1) + 0.125*vertices[23]*(-x[0] + 1)*(x[1] + 1)
diff -r 5b37f431356035e27bca10a03dd4a7b78df0c980 -r 3291782eb19463320b8fdccba398fe66e98d87a1 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -24,7 +24,8 @@
from sympy import \
symarray, \
symbols, \
- diff
+ diff, \
+ ccode
import re
import yaml
@@ -123,14 +124,14 @@
return 'vertices[%d]' % n
def _get_function_line(self, i):
- line = str(self.f[i])
+ line = ccode(self.f[i])
for j in range(self.num_dim):
line = re.sub(r'x_%d' % j, 'x[%d]' % j, line)
line = re.sub(r'(vertices_._.)', self._replace_func, line)
return ''' fx[%d] = %s \n''' % (i, line)
def _get_jacobian_line(self, i, j):
- line = str(self.J[i, j])
+ line = ccode(self.J[i, j])
for k in range(self.num_dim):
line = re.sub(r'x_%d' % k, 'x[%d]' % k, line)
line = re.sub(r'(vertices_._.)', self._replace_func, line)
diff -r 5b37f431356035e27bca10a03dd4a7b78df0c980 -r 3291782eb19463320b8fdccba398fe66e98d87a1 yt/utilities/mesh_types.yaml
--- a/yt/utilities/mesh_types.yaml
+++ b/yt/utilities/mesh_types.yaml
@@ -4,14 +4,14 @@
num_vertices: 8
num_mapped_coords: 3
shape_functions: |
- [(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,
- (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,
- (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.,
+ (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.,
+ (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.,
+ (1 - x[0])*(1 + x[1])*(1 + x[2])/8.]
Quad4:
mesh_type: Q1
@@ -19,10 +19,10 @@
num_vertices: 4
num_mapped_coords: 2
shape_functions: |
- [(1 - x[0])*(1 - x[1])/4,
- (1 + x[0])*(1 - x[1])/4,
- (1 + x[0])*(1 + x[1])/4,
- (1 - x[0])*(1 + x[1])/4]
+ [(1 - x[0])*(1 - x[1])/4.,
+ (1 + x[0])*(1 - x[1])/4.,
+ (1 + x[0])*(1 + x[1])/4.,
+ (1 - x[0])*(1 + x[1])/4.]
Wedge6:
mesh_type: W1
@@ -30,9 +30,9 @@
num_vertices: 6
num_mapped_coords: 3
shape_functions: |
- [(1 - x[0] - x[1]) * (1 - x[2]) / 2,
- x[0] * (1 - x[2]) / 2,
- 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
+ [(1 - x[0] - x[1]) * (1 - x[2]) / 2.,
+ x[0] * (1 - x[2]) / 2.,
+ 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
https://bitbucket.org/yt_analysis/yt/commits/bfbfac6a208e/
Changeset: bfbfac6a208e
Branch: yt
User: atmyers
Date: 2016-08-03 21:04:38+00:00
Summary: Changing the way 2D Jacobians are handled to make them play more nicely with codegen.
Affected #: 2 files
diff -r 3291782eb19463320b8fdccba398fe66e98d87a1 -r bfbfac6a208ef40c597ebdd067a9b891ae3bb180 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -117,10 +117,11 @@
#
# outputs:
#
-# A - A flattened array storing the Jacobian matrix
-# The order of this array is [J11, J12, J21, J22]
+# rcol - the first column of the jacobian
+# scol - the second column of the jacobian
#
-ctypedef void (*jac_type2D)(double* A,
+ctypedef void (*jac_type2D)(double* rcol,
+ double* scol,
double* x,
double* vertices,
double* phys_x) nogil
diff -r 3291782eb19463320b8fdccba398fe66e98d87a1 -r bfbfac6a208ef40c597ebdd067a9b891ae3bb180 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -731,11 +731,11 @@
# begin Newton iteration
while (err > self.tolerance and iterations < self.max_iter):
- self.jac(A, x, vertices, physical_x)
+ 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[1]*f[1]) / d
- x[1] -= (-A[2]*f[0] + A[0]*f[1]) / d
+ 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)
err = maxnorm(f, 2)
https://bitbucket.org/yt_analysis/yt/commits/a94510673133/
Changeset: a94510673133
Branch: yt
User: atmyers
Date: 2016-08-03 21:06:17+00:00
Summary: Changing the code generation to use MatrixSymbol to avoid string pattern matching.
Affected #: 3 files
diff -r bfbfac6a208ef40c597ebdd067a9b891ae3bb180 -r a94510673133ed88a27fa16b910a61d3b247f44b yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -1,10 +1,12 @@
-cdef void Q1Function2D(double* fx,
+cdef void Q1Function3D(double* fx,
double* x,
double* vertices,
double* phys_x) nogil
-cdef void Q1Jacobian2D(double* A,
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
double* x,
double* vertices,
double* phys_x) nogil
@@ -24,15 +26,14 @@
double* phys_x) nogil
-cdef void Q1Function3D(double* fx,
+cdef void Q1Function2D(double* fx,
double* x,
double* vertices,
double* phys_x) nogil
-cdef void Q1Jacobian3D(double* rcol,
+cdef void Q1Jacobian2D(double* rcol,
double* scol,
- double* tcol,
double* x,
double* vertices,
double* phys_x) nogil
diff -r bfbfac6a208ef40c597ebdd067a9b891ae3bb180 -r a94510673133ed88a27fa16b910a61d3b247f44b yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -8,25 +8,33 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void Q1Function2D(double* fx,
+cdef void Q1Function3D(double* fx,
double* x,
double* vertices,
double* phys_x) nogil:
- fx[0] = -phys_x[0] + 0.25*vertices[0]*(-x[0] + 1)*(-x[1] + 1) + 0.25*vertices[2]*(x[0] + 1)*(-x[1] + 1) + 0.25*vertices[4]*(x[0] + 1)*(x[1] + 1) + 0.25*vertices[6]*(-x[0] + 1)*(x[1] + 1)
- fx[1] = -phys_x[1] + 0.25*vertices[1]*(-x[0] + 1)*(-x[1] + 1) + 0.25*vertices[3]*(x[0] + 1)*(-x[1] + 1) + 0.25*vertices[5]*(x[0] + 1)*(x[1] + 1) + 0.25*vertices[7]*(-x[0] + 1)*(x[1] + 1)
+ fx[0] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[21] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[18] - phys_x[0];
+ fx[1] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[22] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[19] - phys_x[1];
+ fx[2] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[23] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[20] - phys_x[2];
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void Q1Jacobian2D(double* A,
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
double* x,
double* vertices,
double* phys_x) nogil:
- A[0] = -0.25*vertices[0]*(-x[1] + 1) + 0.25*vertices[2]*(-x[1] + 1) + 0.25*vertices[4]*(x[1] + 1) - 0.25*vertices[6]*(x[1] + 1)
- A[1] = -0.25*vertices[0]*(-x[0] + 1) - 0.25*vertices[2]*(x[0] + 1) + 0.25*vertices[4]*(x[0] + 1) + 0.25*vertices[6]*(-x[0] + 1)
- A[2] = -0.25*vertices[1]*(-x[1] + 1) + 0.25*vertices[3]*(-x[1] + 1) + 0.25*vertices[5]*(x[1] + 1) - 0.25*vertices[7]*(x[1] + 1)
- A[3] = -0.25*vertices[1]*(-x[0] + 1) - 0.25*vertices[3]*(x[0] + 1) + 0.25*vertices[5]*(x[0] + 1) + 0.25*vertices[7]*(-x[0] + 1)
+ rcol[0] = -0.125*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[1])*(1 - x[2])*vertices[3] - 0.125*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[1])*(1 - x[2])*vertices[6] - 0.125*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 + x[1])*(1 + x[2])*vertices[18] - 0.125*(1 + x[1])*(1 + x[2])*vertices[21];
+ scol[0] = -0.125*(1 - x[0])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[2])*vertices[9] - 0.125*(1 - x[0])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[2])*vertices[21] - 0.125*(1 + x[0])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[2])*vertices[6] - 0.125*(1 + x[0])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[2])*vertices[18];
+ tcol[0] = -0.125*(1 - x[0])*(1 - x[1])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*vertices[12] - 0.125*(1 - x[0])*(1 + x[1])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*vertices[21] - 0.125*(1 + x[0])*(1 - x[1])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*vertices[15] - 0.125*(1 + x[0])*(1 + x[1])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*vertices[18];
+ rcol[1] = -0.125*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[1])*(1 - x[2])*vertices[4] - 0.125*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[1])*(1 - x[2])*vertices[7] - 0.125*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 + x[1])*(1 + x[2])*vertices[19] - 0.125*(1 + x[1])*(1 + x[2])*vertices[22];
+ scol[1] = -0.125*(1 - x[0])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[2])*vertices[10] - 0.125*(1 - x[0])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[2])*vertices[22] - 0.125*(1 + x[0])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[2])*vertices[7] - 0.125*(1 + x[0])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[2])*vertices[19];
+ tcol[1] = -0.125*(1 - x[0])*(1 - x[1])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*vertices[13] - 0.125*(1 - x[0])*(1 + x[1])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*vertices[22] - 0.125*(1 + x[0])*(1 - x[1])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*vertices[16] - 0.125*(1 + x[0])*(1 + x[1])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*vertices[19];
+ rcol[2] = -0.125*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[1])*(1 - x[2])*vertices[5] - 0.125*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[1])*(1 - x[2])*vertices[8] - 0.125*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 + x[1])*(1 + x[2])*vertices[20] - 0.125*(1 + x[1])*(1 + x[2])*vertices[23];
+ scol[2] = -0.125*(1 - x[0])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[2])*vertices[11] - 0.125*(1 - x[0])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[2])*vertices[23] - 0.125*(1 + x[0])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[2])*vertices[8] - 0.125*(1 + x[0])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[2])*vertices[20];
+ tcol[2] = -0.125*(1 - x[0])*(1 - x[1])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*vertices[14] - 0.125*(1 - x[0])*(1 + x[1])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*vertices[23] - 0.125*(1 + x[0])*(1 - x[1])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*vertices[17] - 0.125*(1 + x[0])*(1 + x[1])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*vertices[20];
@cython.boundscheck(False)
@@ -36,9 +44,9 @@
double* x,
double* vertices,
double* phys_x) nogil:
- fx[0] = -phys_x[0] + 0.5*vertices[0]*(-x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[3]*x[0]*(-x[2] + 1) + 0.5*vertices[6]*x[1]*(-x[2] + 1) + 0.5*vertices[9]*(x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[12]*x[0]*(x[2] + 1) + 0.5*vertices[15]*x[1]*(x[2] + 1)
- fx[1] = -phys_x[1] + 0.5*vertices[1]*(-x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[4]*x[0]*(-x[2] + 1) + 0.5*vertices[7]*x[1]*(-x[2] + 1) + 0.5*vertices[10]*(x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[13]*x[0]*(x[2] + 1) + 0.5*vertices[16]*x[1]*(x[2] + 1)
- fx[2] = -phys_x[2] + 0.5*vertices[2]*(-x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[5]*x[0]*(-x[2] + 1) + 0.5*vertices[8]*x[1]*(-x[2] + 1) + 0.5*vertices[11]*(x[2] + 1)*(-x[0] - x[1] + 1) + 0.5*vertices[14]*x[0]*(x[2] + 1) + 0.5*vertices[17]*x[1]*(x[2] + 1)
+ 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];
@cython.boundscheck(False)
@@ -50,46 +58,39 @@
double* x,
double* vertices,
double* phys_x) nogil:
- rcol[0] = -0.5*vertices[0]*(-x[2] + 1) + 0.5*vertices[3]*(-x[2] + 1) - 0.5*vertices[9]*(x[2] + 1) + 0.5*vertices[12]*(x[2] + 1)
- scol[0] = -0.5*vertices[0]*(-x[2] + 1) + 0.5*vertices[6]*(-x[2] + 1) - 0.5*vertices[9]*(x[2] + 1) + 0.5*vertices[15]*(x[2] + 1)
- tcol[0] = -0.5*vertices[0]*(-x[0] - x[1] + 1) - 0.5*vertices[3]*x[0] - 0.5*vertices[6]*x[1] + 0.5*vertices[9]*(-x[0] - x[1] + 1) + 0.5*vertices[12]*x[0] + 0.5*vertices[15]*x[1]
- rcol[1] = -0.5*vertices[1]*(-x[2] + 1) + 0.5*vertices[4]*(-x[2] + 1) - 0.5*vertices[10]*(x[2] + 1) + 0.5*vertices[13]*(x[2] + 1)
- scol[1] = -0.5*vertices[1]*(-x[2] + 1) + 0.5*vertices[7]*(-x[2] + 1) - 0.5*vertices[10]*(x[2] + 1) + 0.5*vertices[16]*(x[2] + 1)
- tcol[1] = -0.5*vertices[1]*(-x[0] - x[1] + 1) - 0.5*vertices[4]*x[0] - 0.5*vertices[7]*x[1] + 0.5*vertices[10]*(-x[0] - x[1] + 1) + 0.5*vertices[13]*x[0] + 0.5*vertices[16]*x[1]
- rcol[2] = -0.5*vertices[2]*(-x[2] + 1) + 0.5*vertices[5]*(-x[2] + 1) - 0.5*vertices[11]*(x[2] + 1) + 0.5*vertices[14]*(x[2] + 1)
- scol[2] = -0.5*vertices[2]*(-x[2] + 1) + 0.5*vertices[8]*(-x[2] + 1) - 0.5*vertices[11]*(x[2] + 1) + 0.5*vertices[17]*(x[2] + 1)
- tcol[2] = -0.5*vertices[2]*(-x[0] - x[1] + 1) - 0.5*vertices[5]*x[0] - 0.5*vertices[8]*x[1] + 0.5*vertices[11]*(-x[0] - x[1] + 1) + 0.5*vertices[14]*x[0] + 0.5*vertices[17]*x[1]
+ 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];
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void Q1Function3D(double* fx,
+cdef void Q1Function2D(double* fx,
double* x,
double* vertices,
double* phys_x) nogil:
- fx[0] = -phys_x[0] + 0.125*vertices[0]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[3]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[6]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[9]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[12]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[15]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[18]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1) + 0.125*vertices[21]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)
- fx[1] = -phys_x[1] + 0.125*vertices[1]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[4]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[7]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[10]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[13]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[16]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[19]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1) + 0.125*vertices[22]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)
- fx[2] = -phys_x[2] + 0.125*vertices[2]*(-x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[5]*(x[0] + 1)*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[8]*(x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[11]*(-x[0] + 1)*(x[1] + 1)*(-x[2] + 1) + 0.125*vertices[14]*(-x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[17]*(x[0] + 1)*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[20]*(x[0] + 1)*(x[1] + 1)*(x[2] + 1) + 0.125*vertices[23]*(-x[0] + 1)*(x[1] + 1)*(x[2] + 1)
+ fx[0] = 0.25*(1 - x[0])*(1 - x[1])*vertices[0] + 0.25*(1 - x[0])*(1 + x[1])*vertices[6] + 0.25*(1 + x[0])*(1 - x[1])*vertices[2] + 0.25*(1 + x[0])*(1 + x[1])*vertices[4] - phys_x[0];
+ fx[1] = 0.25*(1 - x[0])*(1 - x[1])*vertices[1] + 0.25*(1 - x[0])*(1 + x[1])*vertices[7] + 0.25*(1 + x[0])*(1 - x[1])*vertices[3] + 0.25*(1 + x[0])*(1 + x[1])*vertices[5] - phys_x[1];
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void Q1Jacobian3D(double* rcol,
+cdef void Q1Jacobian2D(double* rcol,
double* scol,
- double* tcol,
double* x,
double* vertices,
double* phys_x) nogil:
- rcol[0] = -0.125*vertices[0]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[3]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[6]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[9]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[12]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[15]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[18]*(x[1] + 1)*(x[2] + 1) - 0.125*vertices[21]*(x[1] + 1)*(x[2] + 1)
- scol[0] = -0.125*vertices[0]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[3]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[6]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[9]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[12]*(-x[0] + 1)*(x[2] + 1) - 0.125*vertices[15]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[18]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[21]*(-x[0] + 1)*(x[2] + 1)
- tcol[0] = -0.125*vertices[0]*(-x[0] + 1)*(-x[1] + 1) - 0.125*vertices[3]*(x[0] + 1)*(-x[1] + 1) - 0.125*vertices[6]*(x[0] + 1)*(x[1] + 1) - 0.125*vertices[9]*(-x[0] + 1)*(x[1] + 1) + 0.125*vertices[12]*(-x[0] + 1)*(-x[1] + 1) + 0.125*vertices[15]*(x[0] + 1)*(-x[1] + 1) + 0.125*vertices[18]*(x[0] + 1)*(x[1] + 1) + 0.125*vertices[21]*(-x[0] + 1)*(x[1] + 1)
- rcol[1] = -0.125*vertices[1]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[4]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[7]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[10]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[13]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[16]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[19]*(x[1] + 1)*(x[2] + 1) - 0.125*vertices[22]*(x[1] + 1)*(x[2] + 1)
- scol[1] = -0.125*vertices[1]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[4]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[7]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[10]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[13]*(-x[0] + 1)*(x[2] + 1) - 0.125*vertices[16]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[19]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[22]*(-x[0] + 1)*(x[2] + 1)
- tcol[1] = -0.125*vertices[1]*(-x[0] + 1)*(-x[1] + 1) - 0.125*vertices[4]*(x[0] + 1)*(-x[1] + 1) - 0.125*vertices[7]*(x[0] + 1)*(x[1] + 1) - 0.125*vertices[10]*(-x[0] + 1)*(x[1] + 1) + 0.125*vertices[13]*(-x[0] + 1)*(-x[1] + 1) + 0.125*vertices[16]*(x[0] + 1)*(-x[1] + 1) + 0.125*vertices[19]*(x[0] + 1)*(x[1] + 1) + 0.125*vertices[22]*(-x[0] + 1)*(x[1] + 1)
- rcol[2] = -0.125*vertices[2]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[5]*(-x[1] + 1)*(-x[2] + 1) + 0.125*vertices[8]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[11]*(x[1] + 1)*(-x[2] + 1) - 0.125*vertices[14]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[17]*(-x[1] + 1)*(x[2] + 1) + 0.125*vertices[20]*(x[1] + 1)*(x[2] + 1) - 0.125*vertices[23]*(x[1] + 1)*(x[2] + 1)
- scol[2] = -0.125*vertices[2]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[5]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[8]*(x[0] + 1)*(-x[2] + 1) + 0.125*vertices[11]*(-x[0] + 1)*(-x[2] + 1) - 0.125*vertices[14]*(-x[0] + 1)*(x[2] + 1) - 0.125*vertices[17]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[20]*(x[0] + 1)*(x[2] + 1) + 0.125*vertices[23]*(-x[0] + 1)*(x[2] + 1)
- tcol[2] = -0.125*vertices[2]*(-x[0] + 1)*(-x[1] + 1) - 0.125*vertices[5]*(x[0] + 1)*(-x[1] + 1) - 0.125*vertices[8]*(x[0] + 1)*(x[1] + 1) - 0.125*vertices[11]*(-x[0] + 1)*(x[1] + 1) + 0.125*vertices[14]*(-x[0] + 1)*(-x[1] + 1) + 0.125*vertices[17]*(x[0] + 1)*(-x[1] + 1) + 0.125*vertices[20]*(x[0] + 1)*(x[1] + 1) + 0.125*vertices[23]*(-x[0] + 1)*(x[1] + 1)
+ rcol[0] = -0.25*(1 - x[1])*vertices[0] + 0.25*(1 - x[1])*vertices[2] + 0.25*(1 + x[1])*vertices[4] - 0.25*(1 + x[1])*vertices[6];
+ scol[0] = -0.25*(1 - x[0])*vertices[0] + 0.25*(1 - x[0])*vertices[6] - 0.25*(1 + x[0])*vertices[2] + 0.25*(1 + x[0])*vertices[4];
+ rcol[1] = -0.25*(1 - x[1])*vertices[1] + 0.25*(1 - x[1])*vertices[3] + 0.25*(1 + x[1])*vertices[5] - 0.25*(1 + x[1])*vertices[7];
+ 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];
diff -r bfbfac6a208ef40c597ebdd067a9b891ae3bb180 -r a94510673133ed88a27fa16b910a61d3b247f44b yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -25,7 +25,9 @@
symarray, \
symbols, \
diff, \
- ccode
+ ccode, \
+ Matrix, \
+ MatrixSymbol
import re
import yaml
@@ -53,7 +55,8 @@
@cython.wraparound(False)
@cython.cdivision(True) \n''' + jac_signature_3D + ': \n'
-jac_signature_2D = '''cdef void %s(double* A,
+jac_signature_2D = '''cdef void %s(double* rcol,
+ double* scol,
double* x,
double* vertices,
double* phys_x) nogil'''
@@ -85,9 +88,9 @@
self.num_vertices = mesh_data['num_vertices']
self.num_mapped_coords = mesh_data['num_mapped_coords']
- x = symarray('x', self.num_mapped_coords)
+ x = MatrixSymbol('x', self.num_mapped_coords, 1)
self.x = x
- self.N = eval(mesh_data['shape_functions'])
+ self.N = Matrix(eval(mesh_data['shape_functions']))
self._compute_jacobian()
def _compute_jacobian(self):
@@ -95,67 +98,52 @@
assert(self.num_vertices == len(self.N))
assert(self.num_dim == self.num_mapped_coords)
- X = symarray('vertices', (self.num_dim, self.num_vertices))
- physical_position = symbols(['phys_x[%s] ' % d for d in '012'[:self.num_dim]])
+ X = MatrixSymbol("vertices", self.num_vertices, self.num_dim)
+ self.fx = MatrixSymbol("fx", self.num_dim, 1)
+ physical_position = MatrixSymbol('phys_x', self.num_dim, 1)
- self.f = X.dot(self.N) - physical_position
+ self.f = (self.N.T * Matrix(X)).T - physical_position
self.J = symarray('J', (self.num_dim, self.num_dim))
for i in range(self.num_dim):
for j, var in enumerate(self.x):
- self.J[i][j] = diff(self.f[i], var)
+ self.J[i][j] = diff(self.f[i, 0], var)
+
+ self.rcol = MatrixSymbol("rcol", self.num_dim, 1)
+ self.scol = MatrixSymbol("scol", self.num_dim, 1)
+ self.tcol = MatrixSymbol("tcol", self.num_dim, 1)
self.function_name = '%sFunction%dD' % (self.mesh_type, self.num_dim)
self.function_header = fun_def_template % self.function_name
self.function_declaration = fun_dec_template % self.function_name
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_declaration = jac_dec_template_3D % self.jacobian_name
+
elif (self.num_dim == 2):
self.jacobian_header = jac_def_template_2D % self.jacobian_name
self.jacobian_declaration = jac_dec_template_2D % self.jacobian_name
- def _replace_func(self, match):
- s = match.group(0)
- i = int(s[-3])
- j = int(s[-1])
- n = self.num_dim*j + i
- return 'vertices[%d]' % n
-
- def _get_function_line(self, i):
- line = ccode(self.f[i])
- for j in range(self.num_dim):
- line = re.sub(r'x_%d' % j, 'x[%d]' % j, line)
- line = re.sub(r'(vertices_._.)', self._replace_func, line)
- return ''' fx[%d] = %s \n''' % (i, line)
-
- def _get_jacobian_line(self, i, j):
- line = ccode(self.J[i, j])
- for k in range(self.num_dim):
- line = re.sub(r'x_%d' % k, 'x[%d]' % k, line)
- line = re.sub(r'(vertices_._.)', self._replace_func, line)
- if (self.num_dim == 2):
- return ''' A[%d] = %s \n''' % (2*i + j, line)
- else:
- assert(self.num_dim == 3)
- col = 'rst'[j]
- return ''' %scol[%d] = %s \n''' % (col, i, line)
-
def get_interpolator_definition(self):
'''
This returns the function definitions for the given mesh type.
'''
+
function_code = self.function_header
- for i in range(self.num_mapped_coords):
- function_code += self._get_function_line(i)
-
+ 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):
- for j in range(self.num_dim):
- jacobian_code += self._get_jacobian_line(i, j)
+ 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:
+ continue
+ jacobian_code += '\t' + ccode(self.J[i,2], self.tcol[i, 0]) + '\n'
return function_code, jacobian_code
https://bitbucket.org/yt_analysis/yt/commits/2a481a25b186/
Changeset: 2a481a25b186
Branch: yt
User: atmyers
Date: 2016-08-03 21:31:41+00:00
Summary: Add note about where to find the mesh generation code.
Affected #: 3 files
diff -r a94510673133ed88a27fa16b910a61d3b247f44b -r 2a481a25b186db45c479c09188bf73c6d8fdbeb6 yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -1,17 +1,3 @@
-cdef void Q1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil
-
-
-cdef void Q1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil
-
-
cdef void W1Function3D(double* fx,
double* x,
double* vertices,
@@ -39,3 +25,17 @@
double* phys_x) nogil
+cdef void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
diff -r a94510673133ed88a27fa16b910a61d3b247f44b -r 2a481a25b186db45c479c09188bf73c6d8fdbeb6 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -1,5 +1,8 @@
# This file contains auto-generated functions for sampling
-# inside finite element solutions for various mesh types.
+# inside finite element solutions for various mesh types.
+# To see how the code generation works in detail, see
+# yt/utilities/mesh_code_generation.py.
+
cimport cython
@@ -8,38 +11,6 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
-cdef void Q1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- fx[0] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[21] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[18] - phys_x[0];
- fx[1] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[22] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[19] - phys_x[1];
- fx[2] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[23] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[20] - phys_x[2];
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void Q1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- rcol[0] = -0.125*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[1])*(1 - x[2])*vertices[3] - 0.125*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[1])*(1 - x[2])*vertices[6] - 0.125*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 + x[1])*(1 + x[2])*vertices[18] - 0.125*(1 + x[1])*(1 + x[2])*vertices[21];
- scol[0] = -0.125*(1 - x[0])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[2])*vertices[9] - 0.125*(1 - x[0])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[2])*vertices[21] - 0.125*(1 + x[0])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[2])*vertices[6] - 0.125*(1 + x[0])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[2])*vertices[18];
- tcol[0] = -0.125*(1 - x[0])*(1 - x[1])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*vertices[12] - 0.125*(1 - x[0])*(1 + x[1])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*vertices[21] - 0.125*(1 + x[0])*(1 - x[1])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*vertices[15] - 0.125*(1 + x[0])*(1 + x[1])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*vertices[18];
- rcol[1] = -0.125*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[1])*(1 - x[2])*vertices[4] - 0.125*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[1])*(1 - x[2])*vertices[7] - 0.125*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 + x[1])*(1 + x[2])*vertices[19] - 0.125*(1 + x[1])*(1 + x[2])*vertices[22];
- scol[1] = -0.125*(1 - x[0])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[2])*vertices[10] - 0.125*(1 - x[0])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[2])*vertices[22] - 0.125*(1 + x[0])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[2])*vertices[7] - 0.125*(1 + x[0])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[2])*vertices[19];
- tcol[1] = -0.125*(1 - x[0])*(1 - x[1])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*vertices[13] - 0.125*(1 - x[0])*(1 + x[1])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*vertices[22] - 0.125*(1 + x[0])*(1 - x[1])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*vertices[16] - 0.125*(1 + x[0])*(1 + x[1])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*vertices[19];
- rcol[2] = -0.125*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[1])*(1 - x[2])*vertices[5] - 0.125*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[1])*(1 - x[2])*vertices[8] - 0.125*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 + x[1])*(1 + x[2])*vertices[20] - 0.125*(1 + x[1])*(1 + x[2])*vertices[23];
- scol[2] = -0.125*(1 - x[0])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[2])*vertices[11] - 0.125*(1 - x[0])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[2])*vertices[23] - 0.125*(1 + x[0])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[2])*vertices[8] - 0.125*(1 + x[0])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[2])*vertices[20];
- tcol[2] = -0.125*(1 - x[0])*(1 - x[1])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*vertices[14] - 0.125*(1 - x[0])*(1 + x[1])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*vertices[23] - 0.125*(1 + x[0])*(1 - x[1])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*vertices[17] - 0.125*(1 + x[0])*(1 + x[1])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*vertices[20];
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
cdef void W1Function3D(double* fx,
double* x,
double* vertices,
@@ -94,3 +65,35 @@
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 Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[21] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[18] - phys_x[0];
+ fx[1] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[22] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[19] - phys_x[1];
+ fx[2] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[23] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[20] - phys_x[2];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.125*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[1])*(1 - x[2])*vertices[3] - 0.125*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[1])*(1 - x[2])*vertices[6] - 0.125*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 + x[1])*(1 + x[2])*vertices[18] - 0.125*(1 + x[1])*(1 + x[2])*vertices[21];
+ scol[0] = -0.125*(1 - x[0])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[2])*vertices[9] - 0.125*(1 - x[0])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[2])*vertices[21] - 0.125*(1 + x[0])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[2])*vertices[6] - 0.125*(1 + x[0])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[2])*vertices[18];
+ tcol[0] = -0.125*(1 - x[0])*(1 - x[1])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*vertices[12] - 0.125*(1 - x[0])*(1 + x[1])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*vertices[21] - 0.125*(1 + x[0])*(1 - x[1])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*vertices[15] - 0.125*(1 + x[0])*(1 + x[1])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*vertices[18];
+ rcol[1] = -0.125*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[1])*(1 - x[2])*vertices[4] - 0.125*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[1])*(1 - x[2])*vertices[7] - 0.125*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 + x[1])*(1 + x[2])*vertices[19] - 0.125*(1 + x[1])*(1 + x[2])*vertices[22];
+ scol[1] = -0.125*(1 - x[0])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[2])*vertices[10] - 0.125*(1 - x[0])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[2])*vertices[22] - 0.125*(1 + x[0])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[2])*vertices[7] - 0.125*(1 + x[0])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[2])*vertices[19];
+ tcol[1] = -0.125*(1 - x[0])*(1 - x[1])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*vertices[13] - 0.125*(1 - x[0])*(1 + x[1])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*vertices[22] - 0.125*(1 + x[0])*(1 - x[1])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*vertices[16] - 0.125*(1 + x[0])*(1 + x[1])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*vertices[19];
+ rcol[2] = -0.125*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[1])*(1 - x[2])*vertices[5] - 0.125*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[1])*(1 - x[2])*vertices[8] - 0.125*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 + x[1])*(1 + x[2])*vertices[20] - 0.125*(1 + x[1])*(1 + x[2])*vertices[23];
+ scol[2] = -0.125*(1 - x[0])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[2])*vertices[11] - 0.125*(1 - x[0])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[2])*vertices[23] - 0.125*(1 + x[0])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[2])*vertices[8] - 0.125*(1 + x[0])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[2])*vertices[20];
+ tcol[2] = -0.125*(1 - x[0])*(1 - x[1])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*vertices[14] - 0.125*(1 - x[0])*(1 + x[1])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*vertices[23] - 0.125*(1 + x[0])*(1 - x[1])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*vertices[17] - 0.125*(1 + x[0])*(1 + x[1])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*vertices[20];
+
+
diff -r a94510673133ed88a27fa16b910a61d3b247f44b -r 2a481a25b186db45c479c09188bf73c6d8fdbeb6 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -66,7 +66,9 @@
@cython.cdivision(True) \n''' + jac_signature_2D + ': \n'
file_header = "# This file contains auto-generated functions for sampling \n" + \
- "# inside finite element solutions for various mesh types."
+ "# inside finite element solutions for various mesh types. \n" + \
+ "# To see how the code generation works in detail, see \n" + \
+ "# yt/utilities/mesh_code_generation.py. \n"
class MeshCodeGenerator:
@@ -125,7 +127,8 @@
elif (self.num_dim == 2):
self.jacobian_header = jac_def_template_2D % self.jacobian_name
- self.jacobian_declaration = jac_dec_template_2D % self.jacobian_name
+ self.jacobian_declaration = jac_dec_template_2D % self.jacobian_name
+
def get_interpolator_definition(self):
'''
https://bitbucket.org/yt_analysis/yt/commits/dfe567bb489d/
Changeset: dfe567bb489d
Branch: yt
User: atmyers
Date: 2016-08-03 21:38:45+00:00
Summary: Iterate over the mesh types in sorted order, so you don't get spurious diffs in the generated files.
Affected #: 3 files
diff -r 2a481a25b186db45c479c09188bf73c6d8fdbeb6 -r dfe567bb489d0ce1cb67c2b38cea1c06e99c3e57 yt/utilities/lib/autogenerated_element_samplers.pxd
--- a/yt/utilities/lib/autogenerated_element_samplers.pxd
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -1,10 +1,10 @@
-cdef void W1Function3D(double* fx,
+cdef void Q1Function3D(double* fx,
double* x,
double* vertices,
double* phys_x) nogil
-cdef void W1Jacobian3D(double* rcol,
+cdef void Q1Jacobian3D(double* rcol,
double* scol,
double* tcol,
double* x,
@@ -25,13 +25,13 @@
double* phys_x) nogil
-cdef void Q1Function3D(double* fx,
+cdef void W1Function3D(double* fx,
double* x,
double* vertices,
double* phys_x) nogil
-cdef void Q1Jacobian3D(double* rcol,
+cdef void W1Jacobian3D(double* rcol,
double* scol,
double* tcol,
double* x,
diff -r 2a481a25b186db45c479c09188bf73c6d8fdbeb6 -r dfe567bb489d0ce1cb67c2b38cea1c06e99c3e57 yt/utilities/lib/autogenerated_element_samplers.pyx
--- a/yt/utilities/lib/autogenerated_element_samplers.pyx
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -11,6 +11,63 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
+cdef void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[21] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[18] - phys_x[0];
+ fx[1] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[22] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[19] - phys_x[1];
+ fx[2] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[23] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[20] - phys_x[2];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.125*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[1])*(1 - x[2])*vertices[3] - 0.125*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[1])*(1 - x[2])*vertices[6] - 0.125*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 + x[1])*(1 + x[2])*vertices[18] - 0.125*(1 + x[1])*(1 + x[2])*vertices[21];
+ scol[0] = -0.125*(1 - x[0])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[2])*vertices[9] - 0.125*(1 - x[0])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[2])*vertices[21] - 0.125*(1 + x[0])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[2])*vertices[6] - 0.125*(1 + x[0])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[2])*vertices[18];
+ tcol[0] = -0.125*(1 - x[0])*(1 - x[1])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*vertices[12] - 0.125*(1 - x[0])*(1 + x[1])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*vertices[21] - 0.125*(1 + x[0])*(1 - x[1])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*vertices[15] - 0.125*(1 + x[0])*(1 + x[1])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*vertices[18];
+ rcol[1] = -0.125*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[1])*(1 - x[2])*vertices[4] - 0.125*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[1])*(1 - x[2])*vertices[7] - 0.125*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 + x[1])*(1 + x[2])*vertices[19] - 0.125*(1 + x[1])*(1 + x[2])*vertices[22];
+ scol[1] = -0.125*(1 - x[0])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[2])*vertices[10] - 0.125*(1 - x[0])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[2])*vertices[22] - 0.125*(1 + x[0])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[2])*vertices[7] - 0.125*(1 + x[0])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[2])*vertices[19];
+ tcol[1] = -0.125*(1 - x[0])*(1 - x[1])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*vertices[13] - 0.125*(1 - x[0])*(1 + x[1])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*vertices[22] - 0.125*(1 + x[0])*(1 - x[1])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*vertices[16] - 0.125*(1 + x[0])*(1 + x[1])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*vertices[19];
+ rcol[2] = -0.125*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[1])*(1 - x[2])*vertices[5] - 0.125*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[1])*(1 - x[2])*vertices[8] - 0.125*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 + x[1])*(1 + x[2])*vertices[20] - 0.125*(1 + x[1])*(1 + x[2])*vertices[23];
+ scol[2] = -0.125*(1 - x[0])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[2])*vertices[11] - 0.125*(1 - x[0])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[2])*vertices[23] - 0.125*(1 + x[0])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[2])*vertices[8] - 0.125*(1 + x[0])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[2])*vertices[20];
+ tcol[2] = -0.125*(1 - x[0])*(1 - x[1])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*vertices[14] - 0.125*(1 - x[0])*(1 + x[1])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*vertices[23] - 0.125*(1 + x[0])*(1 - x[1])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*vertices[17] - 0.125*(1 + x[0])*(1 + x[1])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*vertices[20];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = 0.25*(1 - x[0])*(1 - x[1])*vertices[0] + 0.25*(1 - x[0])*(1 + x[1])*vertices[6] + 0.25*(1 + x[0])*(1 - x[1])*vertices[2] + 0.25*(1 + x[0])*(1 + x[1])*vertices[4] - phys_x[0];
+ fx[1] = 0.25*(1 - x[0])*(1 - x[1])*vertices[1] + 0.25*(1 - x[0])*(1 + x[1])*vertices[7] + 0.25*(1 + x[0])*(1 - x[1])*vertices[3] + 0.25*(1 + x[0])*(1 + x[1])*vertices[5] - phys_x[1];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.25*(1 - x[1])*vertices[0] + 0.25*(1 - x[1])*vertices[2] + 0.25*(1 + x[1])*vertices[4] - 0.25*(1 + x[1])*vertices[6];
+ scol[0] = -0.25*(1 - x[0])*vertices[0] + 0.25*(1 - x[0])*vertices[6] - 0.25*(1 + x[0])*vertices[2] + 0.25*(1 + x[0])*vertices[4];
+ rcol[1] = -0.25*(1 - x[1])*vertices[1] + 0.25*(1 - x[1])*vertices[3] + 0.25*(1 + x[1])*vertices[5] - 0.25*(1 + x[1])*vertices[7];
+ 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 W1Function3D(double* fx,
double* x,
double* vertices,
@@ -40,60 +97,3 @@
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];
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void Q1Function2D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- fx[0] = 0.25*(1 - x[0])*(1 - x[1])*vertices[0] + 0.25*(1 - x[0])*(1 + x[1])*vertices[6] + 0.25*(1 + x[0])*(1 - x[1])*vertices[2] + 0.25*(1 + x[0])*(1 + x[1])*vertices[4] - phys_x[0];
- fx[1] = 0.25*(1 - x[0])*(1 - x[1])*vertices[1] + 0.25*(1 - x[0])*(1 + x[1])*vertices[7] + 0.25*(1 + x[0])*(1 - x[1])*vertices[3] + 0.25*(1 + x[0])*(1 + x[1])*vertices[5] - phys_x[1];
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void Q1Jacobian2D(double* rcol,
- double* scol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- rcol[0] = -0.25*(1 - x[1])*vertices[0] + 0.25*(1 - x[1])*vertices[2] + 0.25*(1 + x[1])*vertices[4] - 0.25*(1 + x[1])*vertices[6];
- scol[0] = -0.25*(1 - x[0])*vertices[0] + 0.25*(1 - x[0])*vertices[6] - 0.25*(1 + x[0])*vertices[2] + 0.25*(1 + x[0])*vertices[4];
- rcol[1] = -0.25*(1 - x[1])*vertices[1] + 0.25*(1 - x[1])*vertices[3] + 0.25*(1 + x[1])*vertices[5] - 0.25*(1 + x[1])*vertices[7];
- 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 Q1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- fx[0] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[21] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[18] - phys_x[0];
- fx[1] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[22] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[19] - phys_x[1];
- fx[2] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[23] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[20] - phys_x[2];
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef void Q1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- rcol[0] = -0.125*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[1])*(1 - x[2])*vertices[3] - 0.125*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[1])*(1 - x[2])*vertices[6] - 0.125*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 + x[1])*(1 + x[2])*vertices[18] - 0.125*(1 + x[1])*(1 + x[2])*vertices[21];
- scol[0] = -0.125*(1 - x[0])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[2])*vertices[9] - 0.125*(1 - x[0])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[2])*vertices[21] - 0.125*(1 + x[0])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[2])*vertices[6] - 0.125*(1 + x[0])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[2])*vertices[18];
- tcol[0] = -0.125*(1 - x[0])*(1 - x[1])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*vertices[12] - 0.125*(1 - x[0])*(1 + x[1])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*vertices[21] - 0.125*(1 + x[0])*(1 - x[1])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*vertices[15] - 0.125*(1 + x[0])*(1 + x[1])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*vertices[18];
- rcol[1] = -0.125*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[1])*(1 - x[2])*vertices[4] - 0.125*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[1])*(1 - x[2])*vertices[7] - 0.125*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 + x[1])*(1 + x[2])*vertices[19] - 0.125*(1 + x[1])*(1 + x[2])*vertices[22];
- scol[1] = -0.125*(1 - x[0])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[2])*vertices[10] - 0.125*(1 - x[0])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[2])*vertices[22] - 0.125*(1 + x[0])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[2])*vertices[7] - 0.125*(1 + x[0])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[2])*vertices[19];
- tcol[1] = -0.125*(1 - x[0])*(1 - x[1])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*vertices[13] - 0.125*(1 - x[0])*(1 + x[1])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*vertices[22] - 0.125*(1 + x[0])*(1 - x[1])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*vertices[16] - 0.125*(1 + x[0])*(1 + x[1])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*vertices[19];
- rcol[2] = -0.125*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[1])*(1 - x[2])*vertices[5] - 0.125*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[1])*(1 - x[2])*vertices[8] - 0.125*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 + x[1])*(1 + x[2])*vertices[20] - 0.125*(1 + x[1])*(1 + x[2])*vertices[23];
- scol[2] = -0.125*(1 - x[0])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[2])*vertices[11] - 0.125*(1 - x[0])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[2])*vertices[23] - 0.125*(1 + x[0])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[2])*vertices[8] - 0.125*(1 + x[0])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[2])*vertices[20];
- tcol[2] = -0.125*(1 - x[0])*(1 - x[1])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*vertices[14] - 0.125*(1 - x[0])*(1 + x[1])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*vertices[23] - 0.125*(1 + x[0])*(1 - x[1])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*vertices[17] - 0.125*(1 + x[0])*(1 + x[1])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*vertices[20];
-
-
diff -r 2a481a25b186db45c479c09188bf73c6d8fdbeb6 -r dfe567bb489d0ce1cb67c2b38cea1c06e99c3e57 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -174,7 +174,7 @@
pyx_file.write("cimport cython \n \n")
pyx_file.write("\n \n")
- for mesh_data in mesh_types.values():
+ for _, mesh_data in sorted(mesh_types.items()):
codegen = MeshCodeGenerator(mesh_data)
function_code, jacobian_code = codegen.get_interpolator_definition()
https://bitbucket.org/yt_analysis/yt/commits/20b0faf6ffcc/
Changeset: 20b0faf6ffcc
Branch: yt
User: atmyers
Date: 2016-08-03 22:40:54+00:00
Summary: Fixing flake8 style issues.
Affected #: 1 file
diff -r dfe567bb489d0ce1cb67c2b38cea1c06e99c3e57 -r 20b0faf6ffcc44584b99eca0d091b91931976fc4 yt/utilities/mesh_code_generation.py
--- a/yt/utilities/mesh_code_generation.py
+++ b/yt/utilities/mesh_code_generation.py
@@ -23,12 +23,10 @@
from sympy import \
symarray, \
- symbols, \
diff, \
ccode, \
Matrix, \
MatrixSymbol
-import re
import yaml
https://bitbucket.org/yt_analysis/yt/commits/1ba5918b38d6/
Changeset: 1ba5918b38d6
Branch: yt
User: jzuhone
Date: 2016-08-24 18:10:25+00:00
Summary: Merged in atmyers/yt (pull request #2319)
Code generation for element types that require Newton-Raphson iteration
Affected #: 9 files
diff -r 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -23,6 +23,7 @@
yt/geometry/particle_smooth.c
yt/geometry/selection_routines.c
yt/utilities/amr_utils.c
+yt/utilities/lib/autogenerated_element_samplers.c
yt/utilities/kdtree/forthonf2c.h
yt/utilities/libconfig_wrapper.c
yt/utilities/spatial/ckdtree.c
diff -r 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 MANIFEST.in
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -3,6 +3,7 @@
include yt/visualization/mapserver/html/leaflet/*.css
include yt/visualization/mapserver/html/leaflet/*.js
include yt/visualization/mapserver/html/leaflet/images/*.png
+include yt/utilities/mesh_types.yaml
exclude scripts/pr_backport.py
recursive-include yt *.py *.pyx *.pxd *.h README* *.txt LICENSE* *.cu
recursive-include doc *.rst *.txt *.py *.ipynb *.png *.jpg *.css *.html
diff -r 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 setup.py
--- a/setup.py
+++ b/setup.py
@@ -114,6 +114,9 @@
["yt/utilities/spatial/ckdtree.pyx"],
include_dirs=["yt/utilities/lib/"],
libraries=std_libs),
+ Extension("yt.utilities.lib.autogenerated_element_samplers",
+ ["yt/utilities/lib/autogenerated_element_samplers.pyx"],
+ include_dirs=["yt/utilities/lib/"]),
Extension("yt.utilities.lib.bitarray",
["yt/utilities/lib/bitarray.pyx"],
libraries=std_libs),
diff -r 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 yt/utilities/lib/autogenerated_element_samplers.pxd
--- /dev/null
+++ b/yt/utilities/lib/autogenerated_element_samplers.pxd
@@ -0,0 +1,41 @@
+cdef void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Q1Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil
+
+
+cdef void Q1Jacobian2D(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 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 yt/utilities/lib/autogenerated_element_samplers.pyx
--- /dev/null
+++ b/yt/utilities/lib/autogenerated_element_samplers.pyx
@@ -0,0 +1,99 @@
+# This file contains auto-generated functions for sampling
+# inside finite element solutions for various mesh types.
+# To see how the code generation works in detail, see
+# yt/utilities/mesh_code_generation.py.
+
+
+cimport cython
+
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Function3D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[21] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[18] - phys_x[0];
+ fx[1] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[22] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[19] - phys_x[1];
+ fx[2] = 0.125*(1 - x[0])*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*(1 + x[2])*vertices[23] + 0.125*(1 + x[0])*(1 - x[1])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[1])*(1 - x[2])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*(1 + x[2])*vertices[20] - phys_x[2];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Jacobian3D(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.125*(1 - x[1])*(1 - x[2])*vertices[0] + 0.125*(1 - x[1])*(1 - x[2])*vertices[3] - 0.125*(1 - x[1])*(1 + x[2])*vertices[12] + 0.125*(1 - x[1])*(1 + x[2])*vertices[15] + 0.125*(1 + x[1])*(1 - x[2])*vertices[6] - 0.125*(1 + x[1])*(1 - x[2])*vertices[9] + 0.125*(1 + x[1])*(1 + x[2])*vertices[18] - 0.125*(1 + x[1])*(1 + x[2])*vertices[21];
+ scol[0] = -0.125*(1 - x[0])*(1 - x[2])*vertices[0] + 0.125*(1 - x[0])*(1 - x[2])*vertices[9] - 0.125*(1 - x[0])*(1 + x[2])*vertices[12] + 0.125*(1 - x[0])*(1 + x[2])*vertices[21] - 0.125*(1 + x[0])*(1 - x[2])*vertices[3] + 0.125*(1 + x[0])*(1 - x[2])*vertices[6] - 0.125*(1 + x[0])*(1 + x[2])*vertices[15] + 0.125*(1 + x[0])*(1 + x[2])*vertices[18];
+ tcol[0] = -0.125*(1 - x[0])*(1 - x[1])*vertices[0] + 0.125*(1 - x[0])*(1 - x[1])*vertices[12] - 0.125*(1 - x[0])*(1 + x[1])*vertices[9] + 0.125*(1 - x[0])*(1 + x[1])*vertices[21] - 0.125*(1 + x[0])*(1 - x[1])*vertices[3] + 0.125*(1 + x[0])*(1 - x[1])*vertices[15] - 0.125*(1 + x[0])*(1 + x[1])*vertices[6] + 0.125*(1 + x[0])*(1 + x[1])*vertices[18];
+ rcol[1] = -0.125*(1 - x[1])*(1 - x[2])*vertices[1] + 0.125*(1 - x[1])*(1 - x[2])*vertices[4] - 0.125*(1 - x[1])*(1 + x[2])*vertices[13] + 0.125*(1 - x[1])*(1 + x[2])*vertices[16] + 0.125*(1 + x[1])*(1 - x[2])*vertices[7] - 0.125*(1 + x[1])*(1 - x[2])*vertices[10] + 0.125*(1 + x[1])*(1 + x[2])*vertices[19] - 0.125*(1 + x[1])*(1 + x[2])*vertices[22];
+ scol[1] = -0.125*(1 - x[0])*(1 - x[2])*vertices[1] + 0.125*(1 - x[0])*(1 - x[2])*vertices[10] - 0.125*(1 - x[0])*(1 + x[2])*vertices[13] + 0.125*(1 - x[0])*(1 + x[2])*vertices[22] - 0.125*(1 + x[0])*(1 - x[2])*vertices[4] + 0.125*(1 + x[0])*(1 - x[2])*vertices[7] - 0.125*(1 + x[0])*(1 + x[2])*vertices[16] + 0.125*(1 + x[0])*(1 + x[2])*vertices[19];
+ tcol[1] = -0.125*(1 - x[0])*(1 - x[1])*vertices[1] + 0.125*(1 - x[0])*(1 - x[1])*vertices[13] - 0.125*(1 - x[0])*(1 + x[1])*vertices[10] + 0.125*(1 - x[0])*(1 + x[1])*vertices[22] - 0.125*(1 + x[0])*(1 - x[1])*vertices[4] + 0.125*(1 + x[0])*(1 - x[1])*vertices[16] - 0.125*(1 + x[0])*(1 + x[1])*vertices[7] + 0.125*(1 + x[0])*(1 + x[1])*vertices[19];
+ rcol[2] = -0.125*(1 - x[1])*(1 - x[2])*vertices[2] + 0.125*(1 - x[1])*(1 - x[2])*vertices[5] - 0.125*(1 - x[1])*(1 + x[2])*vertices[14] + 0.125*(1 - x[1])*(1 + x[2])*vertices[17] + 0.125*(1 + x[1])*(1 - x[2])*vertices[8] - 0.125*(1 + x[1])*(1 - x[2])*vertices[11] + 0.125*(1 + x[1])*(1 + x[2])*vertices[20] - 0.125*(1 + x[1])*(1 + x[2])*vertices[23];
+ scol[2] = -0.125*(1 - x[0])*(1 - x[2])*vertices[2] + 0.125*(1 - x[0])*(1 - x[2])*vertices[11] - 0.125*(1 - x[0])*(1 + x[2])*vertices[14] + 0.125*(1 - x[0])*(1 + x[2])*vertices[23] - 0.125*(1 + x[0])*(1 - x[2])*vertices[5] + 0.125*(1 + x[0])*(1 - x[2])*vertices[8] - 0.125*(1 + x[0])*(1 + x[2])*vertices[17] + 0.125*(1 + x[0])*(1 + x[2])*vertices[20];
+ tcol[2] = -0.125*(1 - x[0])*(1 - x[1])*vertices[2] + 0.125*(1 - x[0])*(1 - x[1])*vertices[14] - 0.125*(1 - x[0])*(1 + x[1])*vertices[11] + 0.125*(1 - x[0])*(1 + x[1])*vertices[23] - 0.125*(1 + x[0])*(1 - x[1])*vertices[5] + 0.125*(1 + x[0])*(1 - x[1])*vertices[17] - 0.125*(1 + x[0])*(1 + x[1])*vertices[8] + 0.125*(1 + x[0])*(1 + x[1])*vertices[20];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Function2D(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ fx[0] = 0.25*(1 - x[0])*(1 - x[1])*vertices[0] + 0.25*(1 - x[0])*(1 + x[1])*vertices[6] + 0.25*(1 + x[0])*(1 - x[1])*vertices[2] + 0.25*(1 + x[0])*(1 + x[1])*vertices[4] - phys_x[0];
+ fx[1] = 0.25*(1 - x[0])*(1 - x[1])*vertices[1] + 0.25*(1 - x[0])*(1 + x[1])*vertices[7] + 0.25*(1 + x[0])*(1 - x[1])*vertices[3] + 0.25*(1 + x[0])*(1 + x[1])*vertices[5] - phys_x[1];
+
+
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True)
+cdef void Q1Jacobian2D(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil:
+ rcol[0] = -0.25*(1 - x[1])*vertices[0] + 0.25*(1 - x[1])*vertices[2] + 0.25*(1 + x[1])*vertices[4] - 0.25*(1 + x[1])*vertices[6];
+ scol[0] = -0.25*(1 - x[0])*vertices[0] + 0.25*(1 - x[0])*vertices[6] - 0.25*(1 + x[0])*vertices[2] + 0.25*(1 + x[0])*vertices[4];
+ rcol[1] = -0.25*(1 - x[1])*vertices[1] + 0.25*(1 - x[1])*vertices[3] + 0.25*(1 + x[1])*vertices[5] - 0.25*(1 + x[1])*vertices[7];
+ 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 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 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 yt/utilities/lib/element_mappings.pxd
--- a/yt/utilities/lib/element_mappings.pxd
+++ b/yt/utilities/lib/element_mappings.pxd
@@ -117,10 +117,11 @@
#
# outputs:
#
-# A - A flattened array storing the Jacobian matrix
-# The order of this array is [J11, J12, J21, J22]
+# rcol - the first column of the jacobian
+# scol - the second column of the jacobian
#
-ctypedef void (*jac_type2D)(double* A,
+ctypedef void (*jac_type2D)(double* rcol,
+ double* scol,
double* x,
double* vertices,
double* phys_x) nogil
diff -r 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 yt/utilities/lib/element_mappings.pyx
--- a/yt/utilities/lib/element_mappings.pyx
+++ b/yt/utilities/lib/element_mappings.pyx
@@ -19,6 +19,13 @@
cimport cython
import numpy as np
from libc.math cimport fabs
+from yt.utilities.lib.autogenerated_element_samplers cimport \
+ Q1Function3D, \
+ Q1Jacobian3D, \
+ Q1Function2D, \
+ Q1Jacobian2D, \
+ W1Function3D, \
+ W1Jacobian3D
cdef extern from "platform_dep.h":
double fmax(double x, double y) nogil
@@ -392,70 +399,6 @@
return -1
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- cdef double rm, rp, sm, sp, tm, tp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
- tm = 1.0 - x[2]
- tp = 1.0 + x[2]
-
- for i in range(3):
- fx[i] = vertices[0 + i]*rm*sm*tm \
- + vertices[3 + i]*rp*sm*tm \
- + vertices[6 + i]*rp*sp*tm \
- + vertices[9 + i]*rm*sp*tm \
- + vertices[12 + i]*rm*sm*tp \
- + vertices[15 + i]*rp*sm*tp \
- + vertices[18 + i]*rp*sp*tp \
- + vertices[21 + i]*rm*sp*tp \
- - 8.0*phys_x[i]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- cdef double rm, rp, sm, sp, tm, tp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
- tm = 1.0 - x[2]
- tp = 1.0 + x[2]
-
- for i in range(3):
- rcol[i] = -sm*tm*vertices[0 + i] + sm*tm*vertices[3 + i] + \
- sp*tm*vertices[6 + i] - sp*tm*vertices[9 + i] - \
- sm*tp*vertices[12 + i] + sm*tp*vertices[15 + i] + \
- sp*tp*vertices[18 + i] - sp*tp*vertices[21 + i]
- scol[i] = -rm*tm*vertices[0 + i] - rp*tm*vertices[3 + i] + \
- rp*tm*vertices[6 + i] + rm*tm*vertices[9 + i] - \
- rm*tp*vertices[12 + i] - rp*tp*vertices[15 + i] + \
- rp*tp*vertices[18 + i] + rm*tp*vertices[21 + i]
- tcol[i] = -rm*sm*vertices[0 + i] - rp*sm*vertices[3 + i] - \
- rp*sp*vertices[6 + i] - rm*sp*vertices[9 + i] + \
- rm*sm*vertices[12 + i] + rp*sm*vertices[15 + i] + \
- rp*sp*vertices[18 + i] + rm*sp*vertices[21 + i]
-
-
cdef class S2Sampler3D(NonlinearSolveSampler3D):
'''
@@ -745,52 +688,6 @@
return -1
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void W1Function3D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- for i in range(3):
- fx[i] = vertices[0 + i]*(1.0 - x[0] - x[1])*(1.0 - x[2]) \
- + vertices[3 + i]*x[0]*(1.0 - x[2]) \
- + vertices[6 + i]*x[1]*(1.0 - x[2]) \
- + vertices[9 + i]*(1.0 - x[0] - x[1])*(1.0 + x[2]) \
- + vertices[12 + i]*x[0]*(1.0 + x[2]) \
- + vertices[15 + i]*x[1]*(1.0 + x[2]) \
- - 2.0*phys_x[i]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void W1Jacobian3D(double* rcol,
- double* scol,
- double* tcol,
- double* x,
- double* vertices,
- double* phys_x) nogil:
-
- cdef int i
- for i in range(3):
- rcol[i] = (x[2] - 1.0) * vertices[0 + i] \
- - (x[2] - 1.0) * vertices[3 + i] \
- - (x[2] + 1.0) * vertices[9 + i] \
- + (x[2] + 1.0) * vertices[12 + i]
- scol[i] = (x[2] - 1.0) * vertices[0 + i] \
- - (x[2] - 1.0) * vertices[6 + i] \
- - (x[2] + 1.0) * vertices[9 + i] \
- + (x[2] + 1.0) * vertices[15 + i]
- tcol[i] = (x[0] + x[1] - 1.0) * vertices[0 + i] \
- - x[0] * vertices[3 + i] \
- - x[1] * vertices[6 + i] \
- - (x[0] + x[1] - 1.0) * vertices[9 + i] \
- + x[0] * vertices[12 + i] \
- + x[1] * vertices[15 + i]
-
-
cdef class NonlinearSolveSampler2D(ElementSampler):
'''
@@ -834,11 +731,11 @@
# begin Newton iteration
while (err > self.tolerance and iterations < self.max_iter):
- self.jac(A, x, vertices, physical_x)
+ 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[1]*f[1]) / d
- x[1] -= (-A[2]*f[0] + A[0]*f[1]) / d
+ 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)
err = maxnorm(f, 2)
@@ -893,48 +790,6 @@
return 0
return 1
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Jacobian2D(double* A,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef double rm, rp, sm, sp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
-
- A[0] = -sm*vertices[0] + sm*vertices[2] + sp*vertices[4] - sp*vertices[6]
- A[1] = -rm*vertices[0] - rp*vertices[2] + rp*vertices[4] + rm*vertices[6]
- A[2] = -sm*vertices[1] + sm*vertices[3] + sp*vertices[5] - sp*vertices[7]
- A[3] = -rm*vertices[1] - rp*vertices[3] + rp*vertices[5] + rm*vertices[7]
-
-
- at cython.boundscheck(False)
- at cython.wraparound(False)
- at cython.cdivision(True)
-cdef inline void Q1Function2D(double* fx,
- double* x,
- double* vertices,
- double* phys_x) nogil:
- cdef int i
- cdef double rm, rp, sm, sp
-
- rm = 1.0 - x[0]
- rp = 1.0 + x[0]
- sm = 1.0 - x[1]
- sp = 1.0 + x[1]
-
- for i in range(2):
- fx[i] = vertices[0 + i]*rm*sm \
- + vertices[2 + i]*rp*sm \
- + vertices[4 + i]*rp*sp \
- + vertices[6 + i]*rm*sp \
- - 4.0*phys_x[i]
-
@cython.boundscheck(False)
@cython.wraparound(False)
diff -r 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 yt/utilities/mesh_code_generation.py
--- /dev/null
+++ b/yt/utilities/mesh_code_generation.py
@@ -0,0 +1,192 @@
+"""
+This file contains code for automatically generating the functions and jacobians
+used when sampling inside the supported finite element mesh types. The supported
+mesh types are defined in yt/utilities/mesh_types.yaml.
+
+Usage (from the yt/utilities directory):
+
+python mesh_code_generation.py
+
+This will generate the necessary functions and write them to
+yt/utilities/lib/autogenerated_element_samplers.pyx.
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2016, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from sympy import \
+ symarray, \
+ diff, \
+ ccode, \
+ Matrix, \
+ MatrixSymbol
+import yaml
+
+
+# define some templates used below
+fun_signature = '''cdef void %s(double* fx,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil'''
+
+fun_dec_template = fun_signature + ' \n'
+fun_def_template = '''@cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True) \n''' + fun_signature + ': \n'
+
+jac_signature_3D = '''cdef void %s(double* rcol,
+ double* scol,
+ double* tcol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil'''
+
+jac_dec_template_3D = jac_signature_3D + ' \n'
+jac_def_template_3D = '''@cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True) \n''' + jac_signature_3D + ': \n'
+
+jac_signature_2D = '''cdef void %s(double* rcol,
+ double* scol,
+ double* x,
+ double* vertices,
+ double* phys_x) nogil'''
+jac_dec_template_2D = jac_signature_2D + ' \n'
+jac_def_template_2D = '''@cython.boundscheck(False)
+ at cython.wraparound(False)
+ at cython.cdivision(True) \n''' + jac_signature_2D + ': \n'
+
+file_header = "# This file contains auto-generated functions for sampling \n" + \
+ "# inside finite element solutions for various mesh types. \n" + \
+ "# To see how the code generation works in detail, see \n" + \
+ "# yt/utilities/mesh_code_generation.py. \n"
+
+
+class MeshCodeGenerator:
+ '''
+
+ A class for automatically generating the functions and jacobians used for
+ sampling inside finite element calculations.
+
+ '''
+ def __init__(self, mesh_data):
+ '''
+
+ Mesh data should be a dictionary containing information about the type
+ of elements used. See yt/utilities/mesh_types.yaml for more information.
+
+ '''
+ self.mesh_type = mesh_data['mesh_type']
+ self.num_dim = mesh_data['num_dim']
+ self.num_vertices = mesh_data['num_vertices']
+ self.num_mapped_coords = mesh_data['num_mapped_coords']
+
+ x = MatrixSymbol('x', self.num_mapped_coords, 1)
+ self.x = x
+ self.N = Matrix(eval(mesh_data['shape_functions']))
+ self._compute_jacobian()
+
+ def _compute_jacobian(self):
+
+ assert(self.num_vertices == len(self.N))
+ assert(self.num_dim == self.num_mapped_coords)
+
+ X = MatrixSymbol("vertices", self.num_vertices, self.num_dim)
+ self.fx = MatrixSymbol("fx", self.num_dim, 1)
+ physical_position = MatrixSymbol('phys_x', self.num_dim, 1)
+
+ self.f = (self.N.T * Matrix(X)).T - physical_position
+
+ self.J = symarray('J', (self.num_dim, self.num_dim))
+ for i in range(self.num_dim):
+ for j, var in enumerate(self.x):
+ self.J[i][j] = diff(self.f[i, 0], var)
+
+ self.rcol = MatrixSymbol("rcol", self.num_dim, 1)
+ self.scol = MatrixSymbol("scol", self.num_dim, 1)
+ self.tcol = MatrixSymbol("tcol", self.num_dim, 1)
+
+ self.function_name = '%sFunction%dD' % (self.mesh_type, self.num_dim)
+ self.function_header = fun_def_template % self.function_name
+ self.function_declaration = fun_dec_template % self.function_name
+
+ 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_declaration = jac_dec_template_3D % self.jacobian_name
+
+ elif (self.num_dim == 2):
+ self.jacobian_header = jac_def_template_2D % self.jacobian_name
+ self.jacobian_declaration = jac_dec_template_2D % self.jacobian_name
+
+ def get_interpolator_definition(self):
+ '''
+
+ This returns the function definitions for the given mesh type.
+
+ '''
+
+ 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:
+ continue
+ jacobian_code += '\t' + ccode(self.J[i,2], self.tcol[i, 0]) + '\n'
+
+ return function_code, jacobian_code
+
+ def get_interpolator_declaration(self):
+ '''
+
+ This returns the function declarations for the given mesh type.
+
+ '''
+ return self.function_declaration, self.jacobian_declaration
+
+
+if __name__ == "__main__":
+
+ with open('mesh_types.yaml', 'r') as f:
+ lines = f.read()
+
+ mesh_types = yaml.load(lines)
+
+ pxd_file = open("lib/autogenerated_element_samplers.pxd", "w")
+ pyx_file = open("lib/autogenerated_element_samplers.pyx", "w")
+
+ pyx_file.write(file_header)
+ pyx_file.write("\n \n")
+ pyx_file.write("cimport cython \n \n")
+ pyx_file.write("\n \n")
+
+ for _, mesh_data in sorted(mesh_types.items()):
+ codegen = MeshCodeGenerator(mesh_data)
+
+ function_code, jacobian_code = codegen.get_interpolator_definition()
+ function_decl, jacobian_decl = codegen.get_interpolator_declaration()
+
+ pxd_file.write(function_decl)
+ pxd_file.write("\n \n")
+ pxd_file.write(jacobian_decl)
+ pxd_file.write("\n \n")
+
+ pyx_file.write(function_code)
+ pyx_file.write("\n \n")
+ pyx_file.write(jacobian_code)
+ pyx_file.write("\n \n")
+
+ pxd_file.close()
+ pyx_file.close()
diff -r 3215003cc44051d99c126f9bf58e9977e2708332 -r 1ba5918b38d61ac6f9134e4ca918f0bde599b349 yt/utilities/mesh_types.yaml
--- /dev/null
+++ b/yt/utilities/mesh_types.yaml
@@ -0,0 +1,38 @@
+Hex8:
+ mesh_type: Q1
+ num_dim: 3
+ num_vertices: 8
+ num_mapped_coords: 3
+ shape_functions: |
+ [(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.,
+ (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.,
+ (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
+ num_vertices: 4
+ num_mapped_coords: 2
+ shape_functions: |
+ [(1 - x[0])*(1 - x[1])/4.,
+ (1 + x[0])*(1 - x[1])/4.,
+ (1 + x[0])*(1 + x[1])/4.,
+ (1 - x[0])*(1 + x[1])/4.]
+
+Wedge6:
+ mesh_type: W1
+ num_dim: 3
+ num_vertices: 6
+ num_mapped_coords: 3
+ shape_functions: |
+ [(1 - x[0] - x[1]) * (1 - x[2]) / 2.,
+ x[0] * (1 - x[2]) / 2.,
+ 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
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