[yt-svn] commit/yt: jzuhone: Merged in atmyers/yt (pull request #2319)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Aug 24 11:10:53 PDT 2016
1 new commit in yt:
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