From 8ebae44ef0ec61ef9b6e13803ac18b85fcae1ee3 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Tue, 6 Jun 2023 16:27:03 -0600 Subject: [PATCH 1/4] add relax quadratic tri --- .pre-commit-config.yaml | 9 +++ ansys/mapdl/reader/cython/_relaxmidside.pyx | 82 ++++++++++++++------- ansys/mapdl/reader/cython/vtk_support.c | 44 +++++++++-- 3 files changed, 103 insertions(+), 32 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index ea0194cf..d9d11a88 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -32,6 +32,15 @@ repos: # additional_dependencies: [toml] # exclude: "tests/" +- repo: https://github.com/pre-commit/mirrors-clang-format + rev: v14.0.6 + hooks: + - id: clang-format + files: | + (?x)^( + femorph/cython/[^_].*\.c.* + )$ + - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.4.0 hooks: diff --git a/ansys/mapdl/reader/cython/_relaxmidside.pyx b/ansys/mapdl/reader/cython/_relaxmidside.pyx index 88a8ef32..5c6b72f0 100644 --- a/ansys/mapdl/reader/cython/_relaxmidside.pyx +++ b/ansys/mapdl/reader/cython/_relaxmidside.pyx @@ -11,11 +11,12 @@ from libc.stdint cimport int32_t, int64_t # VTK celltypes ctypedef unsigned char uint8 -cdef uint8 VTK_QUADRATIC_QUAD = 23 cdef uint8 VTK_HEXAHEDRON = 12 cdef uint8 VTK_PYRAMID = 14 cdef uint8 VTK_TETRA = 10 cdef uint8 VTK_WEDGE = 13 +cdef uint8 VTK_QUADRATIC_TRIANGLE = 22 +cdef uint8 VTK_QUADRATIC_QUAD = 23 cdef uint8 VTK_QUADRATIC_HEXAHEDRON = 25 cdef uint8 VTK_QUADRATIC_PYRAMID = 27 cdef uint8 VTK_QUADRATIC_TETRA = 24 @@ -25,16 +26,44 @@ cdef uint8 VTK_QUADRATIC_WEDGE = 26 #============================================================================== # Quadratic element relaxation functions #============================================================================== +cdef inline void relax_mid_tri(int64_t [::1] cellarr, int c, double [:, ::1] pts, + double rfac): + """ + Resets the midside nodes of the quadratic quad starting at index c. + + relaxation factor rfac + + The ordering of the three points defining the cell is point ids (0-2,3-5) + where id #3 is the midedge node between points (0,1); id #4 is the midedge + node between points (1,2); and id #5 is the midedge node between points + (2,0). + + """ + cdef int ind0 = cellarr[c + 0] + cdef int ind1 = cellarr[c + 1] + cdef int ind2 = cellarr[c + 2] + cdef int ind3 = cellarr[c + 3] + cdef int ind4 = cellarr[c + 4] + cdef int ind5 = cellarr[c + 5] + + cdef int j + + for j in range(3): + pts[ind3, j] = pts[ind3, j]*(1 - rfac) + (pts[ind0, j] + pts[ind1, j])*0.5*rfac + pts[ind4, j] = pts[ind4, j]*(1 - rfac) + (pts[ind1, j] + pts[ind2, j])*0.5*rfac + pts[ind5, j] = pts[ind5, j]*(1 - rfac) + (pts[ind2, j] + pts[ind0, j])*0.5*rfac + + cdef inline void relax_mid_quad(int64_t [::1] cellarr, int c, double [:, ::1] pts, double rfac): """ Resets the midside nodes of the quadratic quad starting at index c. - + relaxation factor rfac - + midedge nodes between (0,1), (1,2), (1,2), (0,3) - + """ cdef int ind0 = cellarr[c + 0] cdef int ind1 = cellarr[c + 1] @@ -58,14 +87,14 @@ cdef inline void relax_mid_tet(int64_t [::1] cellarr, int c, double [:, ::1] pts double rfac): """ Resets the midside nodes of the tetrahedral starting at index c - + relaxation factor rfac - + midedge nodes between (0,1), (1,2), (2,0), (0,3), (1,3), and (2,3) - + """ - + cdef int ind0 = cellarr[c + 0] cdef int ind1 = cellarr[c + 1] cdef int ind2 = cellarr[c + 2] @@ -91,7 +120,7 @@ cdef inline void relax_mid_tet(int64_t [::1] cellarr, int c, double [:, ::1] pts cdef inline void relax_mid_pyr(int64_t [::1] cellarr, int c, double [:, ::1] pts, double rfac): """ - + 5 (0, 1) 6 (1, 2) 7 (2, 3) @@ -100,7 +129,7 @@ cdef inline void relax_mid_pyr(int64_t [::1] cellarr, int c, double [:, ::1] pts 10(1, 4) 11(2, 4) 12(3, 4) - + """ cdef int ind0 = cellarr[c + 0] @@ -133,7 +162,7 @@ cdef inline void relax_mid_pyr(int64_t [::1] cellarr, int c, double [:, ::1] pts cdef inline void relax_mid_weg(int64_t [::1] cellarr, int c, double [:, ::1] pts, double rfac): """ - + 6 (0,1) 7 (1,2) 8 (2,0) @@ -159,7 +188,7 @@ cdef inline void relax_mid_weg(int64_t [::1] cellarr, int c, double [:, ::1] pts cdef int ind12= cellarr[c + 12] cdef int ind13= cellarr[c + 13] cdef int ind14= cellarr[c + 14] - + cdef int j for j in range(3): @@ -178,7 +207,7 @@ cdef inline void relax_mid_hex(int64_t [::1] cellarr, int c, double [:, ::1] pts double rfac): """ - + 8 (0,1) 9 (1,2) 10 (2,3) @@ -191,9 +220,9 @@ cdef inline void relax_mid_hex(int64_t [::1] cellarr, int c, double [:, ::1] pts 17 (1,5) 18 (2,6) 19 (3,7) - + """ - + cdef int ind0 = cellarr[c + 0] cdef int ind1 = cellarr[c + 1] cdef int ind2 = cellarr[c + 2] @@ -234,7 +263,7 @@ cdef inline void relax_mid_hex(int64_t [::1] cellarr, int c, double [:, ::1] pts def reset_midside(int64_t [::1] cellarr, uint8 [::1] celltypes, int64_t [::1] offset, double [:, ::1] pts): - """Resets positions of midside nodes to directly between edge nodes. + """Resets positions of midside nodes to directly between edge nodes. Parameters ---------- @@ -246,22 +275,23 @@ def reset_midside(int64_t [::1] cellarr, uint8 [::1] celltypes, pts : double [:, ::1] 3D double point array. - + """ - cdef int c, i + cdef int i cdef int ncells = offset.size cdef uint8 celltype for i in range(ncells): celltype = celltypes[i] - c = offset[i] + 1 + if celltype == VTK_QUADRATIC_TRIANGLE: + relax_mid_tri(cellarr, offset[i] + 1, pts, 1) if celltype == VTK_QUADRATIC_QUAD: - relax_mid_quad(cellarr, c, pts, 1) - if celltype == VTK_QUADRATIC_TETRA: - relax_mid_tet(cellarr, c, pts, 1) + relax_mid_quad(cellarr, offset[i] + 1, pts, 1) + elif celltype == VTK_QUADRATIC_TETRA: + relax_mid_tet(cellarr, offset[i] + 1, pts, 1) elif celltype == VTK_QUADRATIC_PYRAMID: - relax_mid_pyr(cellarr, c, pts, 1) + relax_mid_pyr(cellarr, offset[i] + 1, pts, 1) elif celltype == VTK_QUADRATIC_WEDGE: - relax_mid_weg(cellarr, c, pts, 1) - elif celltype == VTK_QUADRATIC_HEXAHEDRON: - relax_mid_hex(cellarr, c, pts, 1) + relax_mid_weg(cellarr, offset[i] + 1, pts, 1) + elif celltype == VTK_QUADRATIC_HEXAHEDRON: + relax_mid_hex(cellarr, offset[i] + 1, pts, 1) diff --git a/ansys/mapdl/reader/cython/vtk_support.c b/ansys/mapdl/reader/cython/vtk_support.c index 8a469c07..cadfd148 100644 --- a/ansys/mapdl/reader/cython/vtk_support.c +++ b/ansys/mapdl/reader/cython/vtk_support.c @@ -321,6 +321,27 @@ void add_tri(bool build_offset, const int *elem, bool quad){ return; } +void add_tri_missing_midside(bool build_offset, const int *elem, const int n_mid){ + add_cell(build_offset, 6, VTK_QUADRATIC_TRIANGLE); + + int i; + int elem_indices[] = {4, 5, 7}; + + // edge nodes + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; + + // add midside nodes and populate any missing midside with -1 + for (i = 0; i < n_mid && i < 3; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[elem_indices[i]]]; + } + for (i=3; i>n_mid; i--){ + vtk_data.cells[vtk_data.loc++] = -1; + } + + return; +} static inline void add_line(bool build_offset, const int *elem, int nnode){ bool is_quad; @@ -511,13 +532,24 @@ int ans_to_vtk(const int nelem, const int *elem, const int *elem_off, add_quad(build_offset, &elem[off], true); } else { - // Any other case. (We should not reach this point) - // Assuming quad. - - // printf(" The type could not be identified. Check vtk_support.c file"); - // printf("Number of elements is %d\n" , nnode_elem); + // Any other case. Possible when missing midside nodes. is_quad = nnode_elem > 5; - add_quad(build_offset, &elem[off], is_quad); + + // Check degenerate triangle + if (elem[off + 2] == elem[off + 3]) { + if (is_quad){ + add_tri_missing_midside(build_offset, &elem[off], 3 - (8 - nnode_elem)); + } else { + add_tri(build_offset, &elem[off], false); + } + } else { + // printf(" The type could not be identified. Check vtk_support.c file"); + // printf("Number of elements is %d\n" , nnode_elem); + + // Assume quad + add_quad(build_offset, &elem[off], is_quad); + } + } break; case 4: // solid From fa12bd2b4b2adcff9d1be79948f8d1a74f412596 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Tue, 6 Jun 2023 16:32:25 -0600 Subject: [PATCH 2/4] pre-commit changes --- .flake8 | 2 +- .github/stale.yml | 2 +- .github/workflows/testing-and-deployment.yml | 2 +- .pre-commit-config.yaml | 5 +- ansys/mapdl/reader/cython/_archive.pyx | 4 +- ansys/mapdl/reader/cython/_cellqual.pyx | 1400 +-- ansys/mapdl/reader/cython/_parsefull.pyx | 24 +- ansys/mapdl/reader/cython/_reader.pyx | 40 +- ansys/mapdl/reader/cython/archive.c | 320 +- ansys/mapdl/reader/cython/binary_reader.cpp | 764 +- ansys/mapdl/reader/cython/parsefull.c | 517 +- ansys/mapdl/reader/cython/parsefull.h | 2 +- ansys/mapdl/reader/cython/reader.c | 271 +- ansys/mapdl/reader/cython/vtk_support.c | 244 +- tests/elements/shell281/MAPDL.inp | 118 +- tests/testfiles/para/para0.txt | 10 +- tests/testfiles/para/para1.txt | 7934 +++++++++--------- tests/testfiles/para/para2.txt | 20 +- 18 files changed, 5812 insertions(+), 5867 deletions(-) diff --git a/.flake8 b/.flake8 index 96244f9b..7a8bf8c9 100644 --- a/.flake8 +++ b/.flake8 @@ -1,6 +1,6 @@ [flake8] exclude = venv, __init__.py, build, doc/source/examples -# To be added after refactoring code to be compliant: E501 +# To be added after refactoring code to be compliant: E501 select = W191, W291, W293, W391, E115, E117, E122, E124, E125, E225, E231, E301, E303, F401, F403 count = True max-complexity = 10 diff --git a/.github/stale.yml b/.github/stale.yml index 9c50a67d..0e8d5d62 100644 --- a/.github/stale.yml +++ b/.github/stale.yml @@ -17,5 +17,5 @@ markComment: > # Comment to post when closing a stale issue. Set to `false` to disable closeComment: > Closed due to inactivity. Feel free to reopen if there is new information - or interest. + or interest. Thank you! \ No newline at end of file diff --git a/.github/workflows/testing-and-deployment.yml b/.github/workflows/testing-and-deployment.yml index 76f9503f..a5796c99 100644 --- a/.github/workflows/testing-and-deployment.yml +++ b/.github/workflows/testing-and-deployment.yml @@ -245,7 +245,7 @@ jobs: # - name: Display MAPDL Logs # if: always() # run: cat log.txt - + Release: if: github.event_name == 'push' && contains(github.ref, 'refs/tags') needs: [doc_build, build, mac_build] diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d9d11a88..020f1986 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -38,7 +38,7 @@ repos: - id: clang-format files: | (?x)^( - femorph/cython/[^_].*\.c.* + ansys/mapdl/reader/cython/[^_].*\.c.* )$ - repo: https://github.com/pre-commit/pre-commit-hooks @@ -46,8 +46,9 @@ repos: hooks: - id: check-merge-conflict - id: debug-statements + - id: trailing-whitespace + exclude: '.*\.(cdb|rst|dat)$' -# this validates our github workflow files - repo: https://github.com/python-jsonschema/check-jsonschema rev: 0.23.1 hooks: diff --git a/ansys/mapdl/reader/cython/_archive.pyx b/ansys/mapdl/reader/cython/_archive.pyx index d130d086..ef3aa69e 100644 --- a/ansys/mapdl/reader/cython/_archive.pyx +++ b/ansys/mapdl/reader/cython/_archive.pyx @@ -49,7 +49,7 @@ def py_write_nblock(filename, const int [::1] node_id, int max_node_id, Double array of node coordinates angles : np.ndarray, optional - + """ # attach the stream to the python file @@ -84,7 +84,7 @@ def py_write_nblock_float(filename, const int [::1] node_id, int max_node_id, Double array of node coordinates angles : np.ndarray, optional - + """ # attach the stream to the python file cdef FILE* cfile = fopen(filename.encode(), mode.encode()) diff --git a/ansys/mapdl/reader/cython/_cellqual.pyx b/ansys/mapdl/reader/cython/_cellqual.pyx index 89d6ea92..1714c635 100644 --- a/ansys/mapdl/reader/cython/_cellqual.pyx +++ b/ansys/mapdl/reader/cython/_cellqual.pyx @@ -110,12 +110,12 @@ cdef inline void tet_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, double rfac) nogil: """ Resets the midside nodes of the tetrahedral starting at index c - + relaxation factor rfac - + midedge nodes between (0,1), (1,2), (2,0), (0,3), (1,3), and (2,3) - + """ cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] @@ -140,7 +140,7 @@ cdef inline void tet_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, cdef inline void pyr_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, double rfac) nogil: """ - + 5 (0, 1) 6 (1, 2) 7 (2, 3) @@ -149,7 +149,7 @@ cdef inline void pyr_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, 10(1, 4) 11(2, 4) 12(3, 4) - + """ cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] @@ -181,7 +181,7 @@ cdef inline void pyr_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, cdef inline void weg_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, double rfac) nogil: """ - + 6 (0,1) 7 (1,2) 8 (2,0) @@ -207,7 +207,7 @@ cdef inline void weg_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, cdef int ind12= cells[c + 12] cdef int ind13= cells[c + 13] cdef int ind14= cells[c + 14] - + cdef int j for j in range(3): @@ -226,7 +226,7 @@ cdef inline void hex_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, double rfac) nogil: """ - + 8 (0,1) 9 (1,2) 10 (2,3) @@ -239,9 +239,9 @@ cdef inline void hex_relax_mid(int64_t [::1] cells, int c, double [:, ::1] pts, 17 (1,5) 18 (2,6) 19 (3,7) - + """ - + cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] cdef int ind2 = cells[c + 2] @@ -373,7 +373,7 @@ tet_dual_ind[3][2] = 3 # ============================================================================ # Linear shape checking functions -# Tetrahedral edges +# Tetrahedral edges cdef int [4][3] tet_edges tet_edges[0][0] = 1 tet_edges[0][1] = 2 @@ -436,7 +436,7 @@ weg_edges[5][0] = 3 weg_edges[5][1] = 4 weg_edges[5][2] = 2 -# populate hex edges +# populate hex edges cdef int [8][3] hex_edges hex_edges[0][0] = 1 hex_edges[0][1] = 3 @@ -487,34 +487,34 @@ cdef inline double repair_weg(double [8][3] temp_fpos, temp_face_cent[0][j] = (cell_points[0][j] + cell_points[1][j] + cell_points[2][j])*div_fac - + # Face 2: (1, 2, 5, 4) for j in range(3): temp_face_cent[1][j] = (cell_points[0][j] + cell_points[1][j] + cell_points[4][j] + cell_points[3][j])*0.25 - + # Face 3: (2, 3, 6, 5) for j in range(3): temp_face_cent[2][j] = (cell_points[1][j] + cell_points[2][j] + cell_points[5][j] + cell_points[4][j])*0.25 - + # Face 4: (3, 1, 4, 6) for j in range(3): temp_face_cent[3][j] = (cell_points[2][j] + cell_points[0][j] + cell_points[3][j] + cell_points[5][j])*0.25 - + # Face 5: (4, 6, 5) for j in range(3): temp_face_cent[4][j] = (cell_points[3][j] + cell_points[5][j] + cell_points[4][j])*div_fac - + #========================================================================== # Dual Element #========================================================================== @@ -524,7 +524,7 @@ cdef inline double repair_weg(double [8][3] temp_fpos, # Compute means for j in range(3): temp_dual_cent[i][j] = (1 - TAU_WEG)*temp_face_cent[weg_dual_ind[i][0]][j] + \ - TAU_WEG*(temp_face_cent[weg_dual_ind[i][1]][j] + + TAU_WEG*(temp_face_cent[weg_dual_ind[i][1]][j] + temp_face_cent[weg_dual_ind[i][2]][j])*0.5 #### Compute normal vector extending from the dual faces #### @@ -532,10 +532,10 @@ cdef inline double repair_weg(double [8][3] temp_fpos, for i in range(6): vector_sub(temp_face_cent[weg_dual_ind[i][2]], temp_face_cent[weg_dual_ind[i][0]], e0) vector_sub(temp_face_cent[weg_dual_ind[i][1]], temp_face_cent[weg_dual_ind[i][0]], e1) - + # Perform cross product to get dual element face normal vector cross(e0, e1, temp_norm[i]) - + # Scale normals normscale = 2.0/sqrt(vector_norm(temp_norm[i])) @@ -587,34 +587,34 @@ cdef inline float repair_weg_float(float [8][3] temp_fpos, temp_face_cent[0][j] = (cell_points[0][j] + cell_points[1][j] + cell_points[2][j])*div_fac - + # Face 2: (1, 2, 5, 4) for j in range(3): temp_face_cent[1][j] = (cell_points[0][j] + cell_points[1][j] + cell_points[4][j] + cell_points[3][j])*0.25 - + # Face 3: (2, 3, 6, 5) for j in range(3): temp_face_cent[2][j] = (cell_points[1][j] + cell_points[2][j] + cell_points[5][j] + cell_points[4][j])*0.25 - + # Face 4: (3, 1, 4, 6) for j in range(3): temp_face_cent[3][j] = (cell_points[2][j] + cell_points[0][j] + cell_points[3][j] + cell_points[5][j])*0.25 - + # Face 5: (4, 6, 5) for j in range(3): temp_face_cent[4][j] = (cell_points[3][j] + cell_points[5][j] + cell_points[4][j])*div_fac - + #========================================================================== # Dual Element #========================================================================== @@ -624,7 +624,7 @@ cdef inline float repair_weg_float(float [8][3] temp_fpos, # Compute means for j in range(3): temp_dual_cent[i][j] = (1 - TAU_WEG)*temp_face_cent[weg_dual_ind[i][0]][j] + \ - TAU_WEG*(temp_face_cent[weg_dual_ind[i][1]][j] + + TAU_WEG*(temp_face_cent[weg_dual_ind[i][1]][j] + temp_face_cent[weg_dual_ind[i][2]][j])*0.5 #### Compute normal vector extending from the dual faces #### @@ -632,10 +632,10 @@ cdef inline float repair_weg_float(float [8][3] temp_fpos, for i in range(6): vector_sub_float(temp_face_cent[weg_dual_ind[i][2]], temp_face_cent[weg_dual_ind[i][0]], e0) vector_sub_float(temp_face_cent[weg_dual_ind[i][1]], temp_face_cent[weg_dual_ind[i][0]], e1) - + # Perform cross product to get dual element face normal vector cross_float(e0, e1, temp_norm[i]) - + # Scale normals normscale = 2.0/sqrt(vector_norm_float(temp_norm[i])) @@ -673,8 +673,8 @@ cdef inline float repair_weg_float(float [8][3] temp_fpos, cdef inline double tet_lin_qual(int64_t [::1] cells, int c, double [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a tetrahedral cell's edge nodes """ - - cdef int indS, ind0, ind1, ind2 + + cdef int indS, ind0, ind1, ind2 cdef double [3] e0 cdef double [3] e1 cdef double [3] e2 @@ -682,34 +682,34 @@ cdef inline double tet_lin_qual(int64_t [::1] cells, int c, double [:, ::1] pts) cdef int i, j cdef double normjac, tnorm cdef double jac = 1.1 - + for i in range(4): indS = cells[c + i] ind0 = cells[c + tet_edges[i][0]] ind1 = cells[c + tet_edges[i][1]] ind2 = cells[c + tet_edges[i][2]] - + for j in range(3): e0[j] = pts[ind0, j] - pts[indS, j] e1[j] = pts[ind1, j] - pts[indS, j] e2[j] = pts[ind2, j] - pts[indS, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + return jac*1.414213562373095 cdef inline float tet_lin_qual_float(int64_t [::1] cells, int c, float [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a tetrahedral cell's edge nodes """ - - cdef int indS, ind0, ind1, ind2 + + cdef int indS, ind0, ind1, ind2 cdef float [3] e0 cdef float [3] e1 cdef float [3] e2 @@ -717,35 +717,35 @@ cdef inline float tet_lin_qual_float(int64_t [::1] cells, int c, cdef int i, j cdef float normjac, tnorm cdef float jac = 1.1 - + for i in range(4): indS = cells[c + i] ind0 = cells[c + tet_edges[i][0]] ind1 = cells[c + tet_edges[i][1]] ind2 = cells[c + tet_edges[i][2]] - + for j in range(3): e0[j] = pts[ind0, j] - pts[indS, j] e1[j] = pts[ind1, j] - pts[indS, j] e2[j] = pts[ind2, j] - pts[indS, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + return jac*1.414213562373095 cdef inline double pyr_lin_qual(int64_t [::1] cells, int c, double [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a pyramid cell's edge nodes """ - - cdef int indS, ind0, ind1, ind2 - + + cdef int indS, ind0, ind1, ind2 + cdef double [3] e0 cdef double [3] e1 cdef double [3] e2 @@ -753,40 +753,40 @@ cdef inline double pyr_lin_qual(int64_t [::1] cells, int c, cdef int i, j cdef double normjac, tnorm cdef double jac = 1.1 - + for i in range(4): indS = cells[c + i] ind0 = cells[c + pyr_edges[i][0]] ind1 = cells[c + pyr_edges[i][1]] ind2 = cells[c + pyr_edges[i][2]] - + for j in range(3): e0[j] = pts[ind0, j] - pts[indS, j] e1[j] = pts[ind1, j] - pts[indS, j] e2[j] = pts[ind2, j] - pts[indS, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + return jac*1.14 cdef inline float pyr_lin_qual_float(int64_t [::1] cells, int c, float [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a pyramid cell's edge nodes """ - cdef int indS, ind0, ind1, ind2 + cdef int indS, ind0, ind1, ind2 cdef float [3] e0 cdef float [3] e1 cdef float [3] e2 cdef int i, j cdef float normjac, tnorm cdef float jac = 1.1 - + for i in range(4): indS = cells[c + i] ind0 = cells[c + pyr_edges[i][0]] @@ -811,9 +811,9 @@ cdef inline float pyr_lin_qual_float(int64_t [::1] cells, int c, cdef inline double weg_lin_qual(int64_t [::1] cells, int c, double [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a wedge cell's edge nodes """ - - cdef int indS, ind0, ind1, ind2 - + + cdef int indS, ind0, ind1, ind2 + cdef double [3] e0 cdef double [3] e1 cdef double [3] e2 @@ -821,20 +821,20 @@ cdef inline double weg_lin_qual(int64_t [::1] cells, int c, cdef int i, j cdef double normjac, tnorm cdef double jac = 1.1 - + for i in range(6): indS = cells[c + i] ind0 = cells[c + weg_edges[i][0]] ind1 = cells[c + weg_edges[i][1]] ind2 = cells[c + weg_edges[i][2]] - + for j in range(3): e0[j] = pts[ind0, j] - pts[indS, j] e1[j] = pts[ind1, j] - pts[indS, j] e2[j] = pts[ind2, j] - pts[indS, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian @@ -847,27 +847,27 @@ cdef inline double weg_lin_qual(int64_t [::1] cells, int c, cdef inline float weg_lin_qual_float(int64_t [::1] cells, int c, float [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a wedge cell's edge nodes """ - cdef int indS, ind0, ind1, ind2 + cdef int indS, ind0, ind1, ind2 cdef float [3] e0 cdef float [3] e1 cdef float [3] e2 cdef int i, j cdef float normjac, tnorm cdef float jac = 1.1 - + for i in range(6): indS = cells[c + i] ind0 = cells[c + weg_edges[i][0]] ind1 = cells[c + weg_edges[i][1]] ind2 = cells[c + weg_edges[i][2]] - + for j in range(3): e0[j] = pts[ind0, j] - pts[indS, j] e1[j] = pts[ind1, j] - pts[indS, j] e2[j] = pts[ind2, j] - pts[indS, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian @@ -879,8 +879,8 @@ cdef inline float weg_lin_qual_float(int64_t [::1] cells, int c, cdef inline double hex_lin_qual(int64_t [::1] cells, int c, double [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a hexahedrals cell's edge nodes """ - cdef int indS, ind0, ind1, ind2 - + cdef int indS, ind0, ind1, ind2 + cdef double [3] e0 cdef double [3] e1 cdef double [3] e2 @@ -888,13 +888,13 @@ cdef inline double hex_lin_qual(int64_t [::1] cells, int c, double [:, ::1] pts) cdef int i, j cdef double normjac, tnorm cdef double jac = 1.1 - + for i in range(8): indS = cells[c + i] ind0 = cells[c + hex_edges[i][0]] ind1 = cells[c + hex_edges[i][1]] ind2 = cells[c + hex_edges[i][2]] - + for j in range(3): e0[j] = pts[ind0, j] - pts[indS, j] e1[j] = pts[ind1, j] - pts[indS, j] @@ -914,8 +914,8 @@ cdef inline double hex_lin_qual(int64_t [::1] cells, int c, double [:, ::1] pts) cdef inline float hex_lin_qual_float(int64_t [::1] cells, int c, float [:, ::1] pts) nogil: """ Returns minimum scaled jacobian of a hexahedrals cell's edge nodes """ - cdef int indS, ind0, ind1, ind2 - + cdef int indS, ind0, ind1, ind2 + cdef float [3] e0 cdef float [3] e1 cdef float [3] e2 @@ -923,13 +923,13 @@ cdef inline float hex_lin_qual_float(int64_t [::1] cells, int c, cdef int i, j cdef float normjac, tnorm cdef float jac = 1.1 - + for i in range(8): indS = cells[c + i] ind0 = cells[c + hex_edges[i][0]] ind1 = cells[c + hex_edges[i][1]] ind2 = cells[c + hex_edges[i][2]] - + for j in range(3): e0[j] = pts[ind0, j] - pts[indS, j] e1[j] = pts[ind1, j] - pts[indS, j] @@ -957,7 +957,7 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, import numpy as np import sympy - + isopts = np.array([[1, 0, 0, 0], # edge node 0 [0, 1, 0, 0], # edge node 1 [0, 0, 1, 0], # edge node 2 @@ -968,9 +968,9 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, [0.5, 0, 0, 0.5], # midside node 7 (between 0 and 3) [0, 0.5, 0, 0.5], # midside node 8 (between 1 and 3) [0, 0, 0.5, 0.5]]) # midside node 9 (between 2 and 3) - + zeta0, zeta1, zeta2, zeta3 = sympy.symbols('zeta0, zeta1, zeta2, zeta3') - + shape_functions = [zeta0*(2*zeta0 - 1), # edge node 0 zeta1*(2*zeta1 - 1), # edge node 1 zeta2*(2*zeta2 - 1), # edge node 2 @@ -992,8 +992,8 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, for i in range(nfun): for j in range(nvar): nprime[j, i] = sympy.diff(shape_functions[i], variables[j]) - - + + # This is the N_prime functions with the isopts subbed into them pre_j = np.empty((nfun*nvar, nfun)) for i in range(nfun): @@ -1003,23 +1003,23 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, (zeta1, isopts[i, 1]), (zeta2, isopts[i, 2]), (zeta3, isopts[i, 3])]) - + subtract first row to get square system for each jacobian - + for example, for edge node 0 the pre_j is: N 0 1 2 3 4 5 6 7 8 9 [ 3 0 0 0 0 0 0 0 0 0] # N/dzeta0 [ 0 -1 0 0 4 0 0 0 0 0] # N/dzeta1 [ 0 0 -1 0 0 0 4 0 0 0] # N/dzeta2 [ 0 0 0 -1 0 0 0 4 0 0] # N/dzeta3 - + # vectors to form jacobian are e0 = 4*Node4 - Node1 - 3*Node0 e1 = 4*Node6 - Node2 - 3*Node0 e2 = 4*Node7 - Node3 - 3*Node0 - - - + + + # example imperfect tetrahedral pts = np.array([[ 8.54211689, 0.3652738 , 3.41738048], [ 8.47148386, 0.36139067, 3.43469913], @@ -1030,15 +1030,15 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, [ 8.52434853, 0.37927608, 3.42764442], [ 8.51203386, 0.38771899, 3.3855228 ], [ 8.47671735, 0.38577743, 3.39418213], - [ 8.4942655 , 0.40172127, 3.39578674]]) - + [ 8.4942655 , 0.40172127, 3.39578674]]) + # Programmatically this was difficult as the first column always had to be # subtracted and it made it difficult to program in consistent edge indices. # Found it best just to print out the indices e0 = np.empty(3) e1 = np.empty(3) e2 = np.empty(3) - + ind0 = 0 ind1 = 1 ind2 = 2 @@ -1049,7 +1049,7 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, ind7 = 7 ind8 = 8 ind9 = 9 - + #============================================================================== # FOLLOWING CODE PRINTED #============================================================================== @@ -1067,8 +1067,8 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, endtxt += ' + {:d}*pts[ind{:d}, j]'.format(k[i], i) elif k[i] > 0: endtxt += ' - {:d}*pts[ind{:d}, j]'.format(k[i], i) - - + + for m in range(1, 4): txt = 'e{:d}[j] ='.format(m - 1) k = pre_j[4*j + m].astype(np.int_) @@ -1082,17 +1082,17 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, elif k[i] < 0: txt += ' - {:d}*pts[ind{:d}, j]'.format(k[i], i) print ' ' + txt + endtxt - print + print """ - + cdef double [3] e0 cdef double [3] e1 cdef double [3] e2 cdef int j cdef double normjac, tnorm - + cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] cdef int ind2 = cells[c + 2] @@ -1111,7 +1111,7 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, e2[j] = 4*pts[ind7, j] - pts[ind3, j] - 3*pts[ind0, j] # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) # Store this as minimum jacobian so far cdef double jac = triple_product(e1, e2, e0)/tnorm @@ -1123,21 +1123,21 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, e2[j] = 4*pts[ind8, j] - pts[ind3, j] - 4*pts[ind4, j] + pts[ind0, j] # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Edge node 2 ######################### for j in range(3): e0[j] = 4*pts[ind5, j] - pts[ind1, j] - 4*pts[ind6, j] + pts[ind0, j] e1[j] = 3*pts[ind2, j] - 4*pts[ind6, j] + pts[ind0, j] e2[j] = 4*pts[ind9, j] - pts[ind3, j] - 4*pts[ind6, j] + pts[ind0, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian @@ -1149,85 +1149,85 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, e0[j] = 4*pts[ind8, j] - pts[ind1, j] - 4*pts[ind7, j] + pts[ind0, j] e1[j] = 4*pts[ind9, j] - pts[ind2, j] - 4*pts[ind7, j] + pts[ind0, j] e2[j] = 3*pts[ind3, j] - 4*pts[ind7, j] + pts[ind0, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 4 ######################### for j in range(3): e0[j] = + pts[ind1, j] + 2*pts[ind4, j] - pts[ind0, j] - 2*pts[ind4, j] e1[j] = - pts[ind2, j] + 2*pts[ind5, j] + 2*pts[ind6, j] - pts[ind0, j] - 2*pts[ind4, j] e2[j] = - pts[ind3, j] + 2*pts[ind7, j] + 2*pts[ind8, j] - pts[ind0, j] - 2*pts[ind4, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 5 ######################### for j in range(3): e0[j] = + pts[ind1, j] + 2*pts[ind5, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind6, j] e1[j] = + pts[ind2, j] + 2*pts[ind5, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind6, j] e2[j] = - pts[ind3, j] + 2*pts[ind8, j] + 2*pts[ind9, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind6, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 6 ######################### for j in range(3): e0[j] = - pts[ind1, j] + 2*pts[ind4, j] + 2*pts[ind5, j] - pts[ind0, j] - 2*pts[ind6, j] e1[j] = + pts[ind2, j] + 2*pts[ind6, j] - pts[ind0, j] - 2*pts[ind6, j] e2[j] = - pts[ind3, j] + 2*pts[ind7, j] + 2*pts[ind9, j] - pts[ind0, j] - 2*pts[ind6, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 7 ######################### for j in range(3): e0[j] = - pts[ind1, j] + 2*pts[ind4, j] + 2*pts[ind8, j] - pts[ind0, j] - 2*pts[ind7, j] e1[j] = - pts[ind2, j] + 2*pts[ind6, j] + 2*pts[ind9, j] - pts[ind0, j] - 2*pts[ind7, j] e2[j] = + pts[ind3, j] + 2*pts[ind7, j] - pts[ind0, j] - 2*pts[ind7, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 8 ######################### for j in range(3): e0[j] = + pts[ind1, j] + 2*pts[ind8, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind7, j] e1[j] = - pts[ind2, j] + 2*pts[ind5, j] + 2*pts[ind9, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind7, j] e2[j] = + pts[ind3, j] + 2*pts[ind8, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind7, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 9 ######################### for j in range(3): e0[j] = - pts[ind1, j] + 2*pts[ind5, j] + 2*pts[ind8, j] + pts[ind0, j] - 2*pts[ind6, j] - 2*pts[ind7, j] @@ -1235,13 +1235,13 @@ cdef inline double tet_quad_qual(int64_t [::1] cells, int c, e2[j] = + pts[ind3, j] + 2*pts[ind9, j] + pts[ind0, j] - 2*pts[ind6, j] - 2*pts[ind7, j] # normalize the determinant of the jacobian - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + # adjust jacobian # return jac*1.414213562373095 return jac @@ -1255,7 +1255,7 @@ cdef inline float tet_quad_qual_float(int64_t [::1] cells, int c, cdef int j cdef float normjac, tnorm - + cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] cdef int ind2 = cells[c + 2] @@ -1274,7 +1274,7 @@ cdef inline float tet_quad_qual_float(int64_t [::1] cells, int c, e2[j] = 4*pts[ind7, j] - pts[ind3, j] - 3*pts[ind0, j] # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) # Store this as minimum jacobian so far cdef float jac = triple_product_float(e1, e2, e0)/tnorm @@ -1286,21 +1286,21 @@ cdef inline float tet_quad_qual_float(int64_t [::1] cells, int c, e2[j] = 4*pts[ind8, j] - pts[ind3, j] - 4*pts[ind4, j] + pts[ind0, j] # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Edge node 2 ######################### for j in range(3): e0[j] = 4*pts[ind5, j] - pts[ind1, j] - 4*pts[ind6, j] + pts[ind0, j] e1[j] = 3*pts[ind2, j] - 4*pts[ind6, j] + pts[ind0, j] e2[j] = 4*pts[ind9, j] - pts[ind3, j] - 4*pts[ind6, j] + pts[ind0, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian @@ -1312,85 +1312,85 @@ cdef inline float tet_quad_qual_float(int64_t [::1] cells, int c, e0[j] = 4*pts[ind8, j] - pts[ind1, j] - 4*pts[ind7, j] + pts[ind0, j] e1[j] = 4*pts[ind9, j] - pts[ind2, j] - 4*pts[ind7, j] + pts[ind0, j] e2[j] = 3*pts[ind3, j] - 4*pts[ind7, j] + pts[ind0, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 4 ######################### for j in range(3): e0[j] = + pts[ind1, j] + 2*pts[ind4, j] - pts[ind0, j] - 2*pts[ind4, j] e1[j] = - pts[ind2, j] + 2*pts[ind5, j] + 2*pts[ind6, j] - pts[ind0, j] - 2*pts[ind4, j] e2[j] = - pts[ind3, j] + 2*pts[ind7, j] + 2*pts[ind8, j] - pts[ind0, j] - 2*pts[ind4, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 5 ######################### for j in range(3): e0[j] = + pts[ind1, j] + 2*pts[ind5, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind6, j] e1[j] = + pts[ind2, j] + 2*pts[ind5, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind6, j] e2[j] = - pts[ind3, j] + 2*pts[ind8, j] + 2*pts[ind9, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind6, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 6 ######################### for j in range(3): e0[j] = - pts[ind1, j] + 2*pts[ind4, j] + 2*pts[ind5, j] - pts[ind0, j] - 2*pts[ind6, j] e1[j] = + pts[ind2, j] + 2*pts[ind6, j] - pts[ind0, j] - 2*pts[ind6, j] e2[j] = - pts[ind3, j] + 2*pts[ind7, j] + 2*pts[ind9, j] - pts[ind0, j] - 2*pts[ind6, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 7 ######################### for j in range(3): e0[j] = - pts[ind1, j] + 2*pts[ind4, j] + 2*pts[ind8, j] - pts[ind0, j] - 2*pts[ind7, j] e1[j] = - pts[ind2, j] + 2*pts[ind6, j] + 2*pts[ind9, j] - pts[ind0, j] - 2*pts[ind7, j] e2[j] = + pts[ind3, j] + 2*pts[ind7, j] - pts[ind0, j] - 2*pts[ind7, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 8 ######################### for j in range(3): e0[j] = + pts[ind1, j] + 2*pts[ind8, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind7, j] e1[j] = - pts[ind2, j] + 2*pts[ind5, j] + 2*pts[ind9, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind7, j] e2[j] = + pts[ind3, j] + 2*pts[ind8, j] + pts[ind0, j] - 2*pts[ind4, j] - 2*pts[ind7, j] - + # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + ######################### Midside node 9 ######################### for j in range(3): e0[j] = - pts[ind1, j] + 2*pts[ind5, j] + 2*pts[ind8, j] + pts[ind0, j] - 2*pts[ind6, j] - 2*pts[ind7, j] @@ -1398,13 +1398,13 @@ cdef inline float tet_quad_qual_float(int64_t [::1] cells, int c, e2[j] = + pts[ind3, j] + 2*pts[ind9, j] + pts[ind0, j] - 2*pts[ind6, j] - 2*pts[ind7, j] # normalize the determinant of the jacobian - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm # Track minimum jacobian if normjac < jac: jac = normjac - + # adjust jacobian return jac*1.414213562373095 @@ -1414,38 +1414,38 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, """ Shape function checking code generated by printing the following code using the quadradic hexahedral indices in hex_quad_qual and replacing them with: - + # Map wedge indices to hexahedral indices # ANSYS [I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y, Z, A, B] pyr_map = [0, 1, 2, 3, 4, 4, 4, 4, 5, 6, 7, 8, 4, 4, 4, 4, 9,10,11,12] for i in range(len(pyr_map)): hex_quad_edges[hex_quad_edges == i] = pyr_map[i] - + for i in range(4): print '# Node {:d}'.format(i) print 'for j in range(3):' print ' e0[j] = 4*pts[ind{:d}, j] - pts[ind{:d}, j] - 3*pts[ind{:d}, j]'.format(hex_quad_edges[i][1][0], hex_quad_edges[i][0][0], i) print ' e1[j] = 4*pts[ind{:d}, j] - pts[ind{:d}, j] - 3*pts[ind{:d}, j]'.format(hex_quad_edges[i][1][1], hex_quad_edges[i][0][1], i) print ' e2[j] = 4*pts[ind{:d}, j] - pts[ind{:d}, j] - 3*pts[ind{:d}, j]'.format(hex_quad_edges[i][1][2], hex_quad_edges[i][0][2], i) - print + print print '# normalize the determinant of the jacobian' - print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' + print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' print 'normjac = triple_product(e1, e2, e0)/tnorm' - print + print print '# Track minimum jacobian' print 'if normjac < jac:' print ' jac = normjac' - print + print print 'print normjac' - - + + # midside nodes along x moving edges for i in [8, 9, 10, 11, 16, 17, 18, 19]: print '# Node {:d}'.format(i) print 'for j in range(3):' print ' e0[j] = pts[ind{:d}, j] - pts[ind{:d}, j]'.format(hex_quad_edges[i][0][0], hex_quad_edges[i][0][1]) print - + print ' e1[j] = 2*pts[ind{:d}, j] + \\'.format(hex_quad_edges[i][1][0]) print ' 2*pts[ind{:d}, j] - \\'.format(hex_quad_edges[i][1][1]) print ' pts[ind{:d}, j] - \\'.format(hex_quad_edges[i][1][2]) @@ -1463,25 +1463,25 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, print ' pts[ind{:d}, j] - \\'.format(hex_quad_edges[i][2][5]) print ' pts[ind{:d}, j] + \\'.format(hex_quad_edges[i][2][6]) print ' pts[ind{:d}, j]'.format(hex_quad_edges[i][2][7]) - print + print print '# normalize the determinant of the jacobian' - print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' + print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' print 'normjac = triple_product(e1, e2, e0)/tnorm' - print + print print '# Track minimum jacobian' print 'if normjac < jac:' print ' jac = normjac' - print - print 'print normjac' - - + print + print 'print normjac' + + """ cdef double [3] e0 cdef double [3] e1 cdef double [3] e2 - cdef int i, j + cdef int i, j cdef double normjac, tnorm cdef int ind0 = cells[c + 0] @@ -1497,64 +1497,64 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, cdef int ind10 = cells[c + 10] cdef int ind11 = cells[c + 11] cdef int ind12 = cells[c + 12] - + # Node 0 for j in range(3): e0[j] = 4*pts[ind5, j] - pts[ind1, j] - 3*pts[ind0, j] e1[j] = 4*pts[ind8, j] - pts[ind3, j] - 3*pts[ind0, j] e2[j] = 4*pts[ind9, j] - pts[ind4, j] - 3*pts[ind0, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) cdef double jac = triple_product(e1, e2, e0)/tnorm - - + + # Node 1 for j in range(3): e0[j] = 4*pts[ind6, j] - pts[ind2, j] - 3*pts[ind1, j] e1[j] = 4*pts[ind5, j] - pts[ind0, j] - 3*pts[ind1, j] e2[j] = 4*pts[ind10, j] - pts[ind4, j] - 3*pts[ind1, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 2 for j in range(3): e0[j] = 4*pts[ind7, j] - pts[ind3, j] - 3*pts[ind2, j] e1[j] = 4*pts[ind6, j] - pts[ind1, j] - 3*pts[ind2, j] e2[j] = 4*pts[ind11, j] - pts[ind4, j] - 3*pts[ind2, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 3 for j in range(3): e0[j] = 4*pts[ind8, j] - pts[ind0, j] - 3*pts[ind3, j] e1[j] = 4*pts[ind7, j] - pts[ind2, j] - 3*pts[ind3, j] e2[j] = 4*pts[ind12, j] - pts[ind4, j] - 3*pts[ind3, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 8 for j in range(3): e0[j] = pts[ind1, j] - pts[ind0, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind8, j] - \ pts[ind0, j] - \ @@ -1563,7 +1563,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind5, j] + \ pts[ind7, j] - + e2[j] = 2*pts[ind9, j] + \ 2*pts[ind10, j] - \ pts[ind0, j] - \ @@ -1572,19 +1572,19 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind5, j] + \ pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 9 for j in range(3): e0[j] = pts[ind2, j] - pts[ind1, j] - + e1[j] = 2*pts[ind5, j] + \ 2*pts[ind7, j] - \ pts[ind0, j] - \ @@ -1593,7 +1593,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind6, j] + \ pts[ind8, j] - + e2[j] = 2*pts[ind10, j] + \ 2*pts[ind11, j] - \ pts[ind1, j] - \ @@ -1602,19 +1602,19 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind6, j] + \ pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 10 for j in range(3): e0[j] = pts[ind3, j] - pts[ind2, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind8, j] - \ pts[ind0, j] - \ @@ -1623,7 +1623,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind7, j] + \ pts[ind5, j] - + e2[j] = 2*pts[ind11, j] + \ 2*pts[ind12, j] - \ pts[ind2, j] - \ @@ -1632,19 +1632,19 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind7, j] + \ pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 11 for j in range(3): e0[j] = pts[ind3, j] - pts[ind0, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind12, j] - \ pts[ind0, j] - \ @@ -1653,7 +1653,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind8, j] + \ pts[ind4, j] - + e2[j] = 2*pts[ind5, j] + \ 2*pts[ind7, j] - \ pts[ind0, j] - \ @@ -1662,19 +1662,19 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind8, j] + \ pts[ind6, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 16 for j in range(3): e0[j] = pts[ind0, j] - pts[ind4, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -1683,7 +1683,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind9, j] + \ pts[ind12, j] - + e2[j] = 2*pts[ind5, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -1692,20 +1692,20 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind9, j] + \ pts[ind10, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 17 for j in range(3): e0[j] = pts[ind1, j] - pts[ind4, j] - + e1[j] = 2*pts[ind5, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -1714,7 +1714,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind10, j] + \ pts[ind9, j] - + e2[j] = 2*pts[ind6, j] + \ 2*pts[ind4, j] - \ pts[ind1, j] - \ @@ -1723,20 +1723,20 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind10, j] + \ pts[ind11, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 18 for j in range(3): e0[j] = pts[ind2, j] - pts[ind4, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind4, j] - \ pts[ind1, j] - \ @@ -1745,7 +1745,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind11, j] + \ pts[ind10, j] - + e2[j] = 2*pts[ind7, j] + \ 2*pts[ind4, j] - \ pts[ind2, j] - \ @@ -1754,19 +1754,19 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind11, j] + \ pts[ind12, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 19 for j in range(3): e0[j] = pts[ind3, j] - pts[ind4, j] - + e1[j] = 2*pts[ind7, j] + \ 2*pts[ind4, j] - \ pts[ind2, j] - \ @@ -1775,7 +1775,7 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind12, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -1784,11 +1784,11 @@ cdef inline double pyr_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind12, j] + \ pts[ind9, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac @@ -1802,7 +1802,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, cdef float [3] e1 cdef float [3] e2 - cdef int i, j + cdef int i, j cdef float normjac, tnorm cdef int ind0 = cells[c + 0] @@ -1818,64 +1818,64 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, cdef int ind10 = cells[c + 10] cdef int ind11 = cells[c + 11] cdef int ind12 = cells[c + 12] - + # Node 0 for j in range(3): e0[j] = 4*pts[ind5, j] - pts[ind1, j] - 3*pts[ind0, j] e1[j] = 4*pts[ind8, j] - pts[ind3, j] - 3*pts[ind0, j] e2[j] = 4*pts[ind9, j] - pts[ind4, j] - 3*pts[ind0, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) cdef float jac = triple_product_float(e1, e2, e0)/tnorm - - + + # Node 1 for j in range(3): e0[j] = 4*pts[ind6, j] - pts[ind2, j] - 3*pts[ind1, j] e1[j] = 4*pts[ind5, j] - pts[ind0, j] - 3*pts[ind1, j] e2[j] = 4*pts[ind10, j] - pts[ind4, j] - 3*pts[ind1, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 2 for j in range(3): e0[j] = 4*pts[ind7, j] - pts[ind3, j] - 3*pts[ind2, j] e1[j] = 4*pts[ind6, j] - pts[ind1, j] - 3*pts[ind2, j] e2[j] = 4*pts[ind11, j] - pts[ind4, j] - 3*pts[ind2, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 3 for j in range(3): e0[j] = 4*pts[ind8, j] - pts[ind0, j] - 3*pts[ind3, j] e1[j] = 4*pts[ind7, j] - pts[ind2, j] - 3*pts[ind3, j] e2[j] = 4*pts[ind12, j] - pts[ind4, j] - 3*pts[ind3, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 8 for j in range(3): e0[j] = pts[ind1, j] - pts[ind0, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind8, j] - \ pts[ind0, j] - \ @@ -1884,7 +1884,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind5, j] + \ pts[ind7, j] - + e2[j] = 2*pts[ind9, j] + \ 2*pts[ind10, j] - \ pts[ind0, j] - \ @@ -1893,19 +1893,19 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind5, j] + \ pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 9 for j in range(3): e0[j] = pts[ind2, j] - pts[ind1, j] - + e1[j] = 2*pts[ind5, j] + \ 2*pts[ind7, j] - \ pts[ind0, j] - \ @@ -1914,7 +1914,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind6, j] + \ pts[ind8, j] - + e2[j] = 2*pts[ind10, j] + \ 2*pts[ind11, j] - \ pts[ind1, j] - \ @@ -1923,19 +1923,19 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind6, j] + \ pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 10 for j in range(3): e0[j] = pts[ind3, j] - pts[ind2, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind8, j] - \ pts[ind0, j] - \ @@ -1944,7 +1944,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind7, j] + \ pts[ind5, j] - + e2[j] = 2*pts[ind11, j] + \ 2*pts[ind12, j] - \ pts[ind2, j] - \ @@ -1953,19 +1953,19 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind7, j] + \ pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 11 for j in range(3): e0[j] = pts[ind3, j] - pts[ind0, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind12, j] - \ pts[ind0, j] - \ @@ -1974,7 +1974,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind8, j] + \ pts[ind4, j] - + e2[j] = 2*pts[ind5, j] + \ 2*pts[ind7, j] - \ pts[ind0, j] - \ @@ -1983,19 +1983,19 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind8, j] + \ pts[ind6, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 16 for j in range(3): e0[j] = pts[ind0, j] - pts[ind4, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -2004,7 +2004,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind9, j] + \ pts[ind12, j] - + e2[j] = 2*pts[ind5, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -2013,20 +2013,20 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind9, j] + \ pts[ind10, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 17 for j in range(3): e0[j] = pts[ind1, j] - pts[ind4, j] - + e1[j] = 2*pts[ind5, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -2035,7 +2035,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind10, j] + \ pts[ind9, j] - + e2[j] = 2*pts[ind6, j] + \ 2*pts[ind4, j] - \ pts[ind1, j] - \ @@ -2044,20 +2044,20 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind10, j] + \ pts[ind11, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 18 for j in range(3): e0[j] = pts[ind2, j] - pts[ind4, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind4, j] - \ pts[ind1, j] - \ @@ -2066,7 +2066,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind11, j] + \ pts[ind10, j] - + e2[j] = 2*pts[ind7, j] + \ 2*pts[ind4, j] - \ pts[ind2, j] - \ @@ -2075,19 +2075,19 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind11, j] + \ pts[ind12, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 19 for j in range(3): e0[j] = pts[ind3, j] - pts[ind4, j] - + e1[j] = 2*pts[ind7, j] + \ 2*pts[ind4, j] - \ pts[ind2, j] - \ @@ -2096,7 +2096,7 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind12, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind4, j] - \ pts[ind0, j] - \ @@ -2105,11 +2105,11 @@ cdef inline float pyr_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind12, j] + \ pts[ind9, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac @@ -2126,9 +2126,9 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, cdef double [3] e1 cdef double [3] e2 - cdef int i, j + cdef int i, j cdef double normjac, tnorm - + cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] cdef int ind2 = cells[c + 2] @@ -2144,14 +2144,14 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, cdef int ind12 = cells[c + 12] cdef int ind13 = cells[c + 13] cdef int ind14 = cells[c + 14] - + # Node 0 for j in range(3): e0[j] = 4*pts[ind6, j] - pts[ind1, j] - 3*pts[ind0, j] e1[j] = 4*pts[ind8, j] - pts[ind2, j] - 3*pts[ind0, j] e2[j] = 4*pts[ind12, j] - pts[ind3, j] - 3*pts[ind0, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) # Store current minimum jacobian cdef double jac = -triple_product(e1, e2, e0)/tnorm @@ -2162,73 +2162,73 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, e1[j] = 4*pts[ind6, j] - pts[ind0, j] - 3*pts[ind1, j] e2[j] = 4*pts[ind13, j] - pts[ind4, j] - 3*pts[ind1, j] - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 2 # Since node 3 doesn't exist, special treatment for node 2 for j in range(3): e0[j] = 4*pts[ind8, j] - pts[ind0, j] - 3*pts[ind2, j] e1[j] = 4*pts[ind7, j] - pts[ind1, j] - 3*pts[ind2, j] e2[j] = 4*pts[ind14, j] - pts[ind5, j] - 3*pts[ind2, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 4 for j in range(3): e0[j] = 4*pts[ind11, j] - pts[ind5, j] - 3*pts[ind3, j] e1[j] = 4*pts[ind9, j] - pts[ind4, j] - 3*pts[ind3, j] e2[j] = 4*pts[ind12, j] - pts[ind0, j] - 3*pts[ind3, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 5 for j in range(3): e0[j] = 4*pts[ind9, j] - pts[ind3, j] - 3*pts[ind4, j] e1[j] = 4*pts[ind10, j] - pts[ind5, j] - 3*pts[ind4, j] e2[j] = 4*pts[ind13, j] - pts[ind1, j] - 3*pts[ind4, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 6 # Since node 3 doesn't exist, special treatment for node 2 for j in range(3): e0[j] = 4*pts[ind10, j] - pts[ind4, j] - 3*pts[ind5, j] e1[j] = 4*pts[ind11, j] - pts[ind3, j] - 3*pts[ind5, j] e2[j] = 4*pts[ind14, j] - pts[ind2, j] - 3*pts[ind5, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 8 for j in range(3): e0[j] = pts[ind1, j] - pts[ind0, j] - + e1[j] = 2*pts[ind7, j] + \ 2*pts[ind8, j] - \ pts[ind0, j] - \ @@ -2237,7 +2237,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind2, j] - \ pts[ind6, j] + \ pts[ind2, j] - + e2[j] = 2*pts[ind12, j] + \ 2*pts[ind13, j] - \ pts[ind0, j] - \ @@ -2246,18 +2246,18 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind6, j] + \ pts[ind9, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 9 for j in range(3): e0[j] = pts[ind2, j] - pts[ind1, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind2, j] - \ pts[ind0, j] - \ @@ -2266,7 +2266,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind2, j] - \ pts[ind7, j] + \ pts[ind8, j] - + e2[j] = 2*pts[ind13, j] + \ 2*pts[ind14, j] - \ pts[ind1, j] - \ @@ -2275,18 +2275,18 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind7, j] + \ pts[ind10, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 11 for j in range(3): e0[j] = pts[ind2, j] - pts[ind0, j] - + e1[j] = 2*pts[ind12, j] + \ 2*pts[ind14, j] - \ pts[ind0, j] - \ @@ -2295,7 +2295,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind8, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind6, j] + \ 2*pts[ind2, j] - \ pts[ind0, j] - \ @@ -2304,18 +2304,18 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind2, j] - \ pts[ind8, j] + \ pts[ind7, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 12 for j in range(3): e0[j] = pts[ind3, j] - pts[ind4, j] - + e1[j] = 2*pts[ind10, j] + \ 2*pts[ind11, j] - \ pts[ind3, j] - \ @@ -2324,7 +2324,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind9, j] + \ pts[ind5, j] - + e2[j] = 2*pts[ind12, j] + \ 2*pts[ind13, j] - \ pts[ind0, j] - \ @@ -2333,18 +2333,18 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind9, j] + \ pts[ind6, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 13 for j in range(3): e0[j] = pts[ind4, j] - pts[ind5, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind5, j] - \ pts[ind3, j] - \ @@ -2353,7 +2353,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind10, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind13, j] + \ 2*pts[ind14, j] - \ pts[ind1, j] - \ @@ -2362,18 +2362,18 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind10, j] + \ pts[ind7, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 15 for j in range(3): e0[j] = pts[ind5, j] - pts[ind3, j] - + e1[j] = 2*pts[ind5, j] + \ 2*pts[ind9, j] - \ pts[ind3, j] - \ @@ -2382,7 +2382,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind11, j] + \ pts[ind10, j] - + e2[j] = 2*pts[ind12, j] + \ 2*pts[ind14, j] - \ pts[ind0, j] - \ @@ -2391,18 +2391,18 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind11, j] + \ pts[ind8, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 16 for j in range(3): e0[j] = pts[ind0, j] - pts[ind3, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -2411,7 +2411,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind12, j] + \ pts[ind14, j] - + e2[j] = 2*pts[ind6, j] + \ 2*pts[ind9, j] - \ pts[ind0, j] - \ @@ -2420,18 +2420,18 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind12, j] + \ pts[ind13, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 17 for j in range(3): e0[j] = pts[ind1, j] - pts[ind4, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind9, j] - \ pts[ind0, j] - \ @@ -2440,7 +2440,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind13, j] + \ pts[ind12, j] - + e2[j] = 2*pts[ind7, j] + \ 2*pts[ind10, j] - \ pts[ind1, j] - \ @@ -2449,19 +2449,19 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind13, j] + \ pts[ind14, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Special treatment for midside node 14 (combine responses from 18 and 19 on hex) # Node 17 for j in range(3): e0[j] = pts[ind2, j] - pts[ind5, j] - + e1[j] = 2*pts[ind7, j] + \ 2*pts[ind10, j] - \ pts[ind1, j] - \ @@ -2470,7 +2470,7 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind14, j] + \ pts[ind13, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -2479,14 +2479,14 @@ cdef inline double weg_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind14, j] + \ pts[ind12, j] - - tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) + + tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = -triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + return jac*1.154 @@ -2499,9 +2499,9 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, cdef float [3] e1 cdef float [3] e2 - cdef int i, j + cdef int i, j cdef float normjac, tnorm - + cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] cdef int ind2 = cells[c + 2] @@ -2517,14 +2517,14 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, cdef int ind12 = cells[c + 12] cdef int ind13 = cells[c + 13] cdef int ind14 = cells[c + 14] - + # Node 0 for j in range(3): e0[j] = 4*pts[ind6, j] - pts[ind1, j] - 3*pts[ind0, j] e1[j] = 4*pts[ind8, j] - pts[ind2, j] - 3*pts[ind0, j] e2[j] = 4*pts[ind12, j] - pts[ind3, j] - 3*pts[ind0, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) # Store current minimum jacobian cdef float jac = -triple_product_float(e1, e2, e0)/tnorm @@ -2535,73 +2535,73 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, e1[j] = 4*pts[ind6, j] - pts[ind0, j] - 3*pts[ind1, j] e2[j] = 4*pts[ind13, j] - pts[ind4, j] - 3*pts[ind1, j] - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 2 # Since node 3 doesn't exist, special treatment for node 2 for j in range(3): e0[j] = 4*pts[ind8, j] - pts[ind0, j] - 3*pts[ind2, j] e1[j] = 4*pts[ind7, j] - pts[ind1, j] - 3*pts[ind2, j] e2[j] = 4*pts[ind14, j] - pts[ind5, j] - 3*pts[ind2, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 4 for j in range(3): e0[j] = 4*pts[ind11, j] - pts[ind5, j] - 3*pts[ind3, j] e1[j] = 4*pts[ind9, j] - pts[ind4, j] - 3*pts[ind3, j] e2[j] = 4*pts[ind12, j] - pts[ind0, j] - 3*pts[ind3, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - - + + # Node 5 for j in range(3): e0[j] = 4*pts[ind9, j] - pts[ind3, j] - 3*pts[ind4, j] e1[j] = 4*pts[ind10, j] - pts[ind5, j] - 3*pts[ind4, j] e2[j] = 4*pts[ind13, j] - pts[ind1, j] - 3*pts[ind4, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 6 # Since node 3 doesn't exist, special treatment for node 2 for j in range(3): e0[j] = 4*pts[ind10, j] - pts[ind4, j] - 3*pts[ind5, j] e1[j] = 4*pts[ind11, j] - pts[ind3, j] - 3*pts[ind5, j] e2[j] = 4*pts[ind14, j] - pts[ind2, j] - 3*pts[ind5, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 8 for j in range(3): e0[j] = pts[ind1, j] - pts[ind0, j] - + e1[j] = 2*pts[ind7, j] + \ 2*pts[ind8, j] - \ pts[ind0, j] - \ @@ -2610,7 +2610,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind2, j] - \ pts[ind6, j] + \ pts[ind2, j] - + e2[j] = 2*pts[ind12, j] + \ 2*pts[ind13, j] - \ pts[ind0, j] - \ @@ -2619,18 +2619,18 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind6, j] + \ pts[ind9, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 9 for j in range(3): e0[j] = pts[ind2, j] - pts[ind1, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind2, j] - \ pts[ind0, j] - \ @@ -2639,7 +2639,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind2, j] - \ pts[ind7, j] + \ pts[ind8, j] - + e2[j] = 2*pts[ind13, j] + \ 2*pts[ind14, j] - \ pts[ind1, j] - \ @@ -2648,18 +2648,18 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind7, j] + \ pts[ind10, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 11 for j in range(3): e0[j] = pts[ind2, j] - pts[ind0, j] - + e1[j] = 2*pts[ind12, j] + \ 2*pts[ind14, j] - \ pts[ind0, j] - \ @@ -2668,7 +2668,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind8, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind6, j] + \ 2*pts[ind2, j] - \ pts[ind0, j] - \ @@ -2677,18 +2677,18 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind2, j] - \ pts[ind8, j] + \ pts[ind7, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 12 for j in range(3): e0[j] = pts[ind3, j] - pts[ind4, j] - + e1[j] = 2*pts[ind10, j] + \ 2*pts[ind11, j] - \ pts[ind3, j] - \ @@ -2697,7 +2697,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind9, j] + \ pts[ind5, j] - + e2[j] = 2*pts[ind12, j] + \ 2*pts[ind13, j] - \ pts[ind0, j] - \ @@ -2706,18 +2706,18 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind9, j] + \ pts[ind6, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 13 for j in range(3): e0[j] = pts[ind4, j] - pts[ind5, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind5, j] - \ pts[ind3, j] - \ @@ -2726,7 +2726,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind10, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind13, j] + \ 2*pts[ind14, j] - \ pts[ind1, j] - \ @@ -2735,18 +2735,18 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind10, j] + \ pts[ind7, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 15 for j in range(3): e0[j] = pts[ind5, j] - pts[ind3, j] - + e1[j] = 2*pts[ind5, j] + \ 2*pts[ind9, j] - \ pts[ind3, j] - \ @@ -2755,7 +2755,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind11, j] + \ pts[ind10, j] - + e2[j] = 2*pts[ind12, j] + \ 2*pts[ind14, j] - \ pts[ind0, j] - \ @@ -2764,18 +2764,18 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind11, j] + \ pts[ind8, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 16 for j in range(3): e0[j] = pts[ind0, j] - pts[ind3, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -2784,7 +2784,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind12, j] + \ pts[ind14, j] - + e2[j] = 2*pts[ind6, j] + \ 2*pts[ind9, j] - \ pts[ind0, j] - \ @@ -2793,18 +2793,18 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind12, j] + \ pts[ind13, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 17 for j in range(3): e0[j] = pts[ind1, j] - pts[ind4, j] - + e1[j] = 2*pts[ind6, j] + \ 2*pts[ind9, j] - \ pts[ind0, j] - \ @@ -2813,7 +2813,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind4, j] - \ pts[ind13, j] + \ pts[ind12, j] - + e2[j] = 2*pts[ind7, j] + \ 2*pts[ind10, j] - \ pts[ind1, j] - \ @@ -2822,19 +2822,19 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind13, j] + \ pts[ind14, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Special treatment for midside node 14 (combine responses from 18 and 19 on hex) # Node 17 for j in range(3): e0[j] = pts[ind2, j] - pts[ind5, j] - + e1[j] = 2*pts[ind7, j] + \ 2*pts[ind10, j] - \ pts[ind1, j] - \ @@ -2843,7 +2843,7 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind14, j] + \ pts[ind13, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -2852,14 +2852,14 @@ cdef inline float weg_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind14, j] + \ pts[ind12, j] - - tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) + + tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = -triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + return jac*1.154 @@ -2868,11 +2868,11 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, """ Quadradic shape functions derived and then printed in pre-index form using the following code: - + import sympy import numpy as np - # idea node locations for quadradic hexahedral + # idea node locations for quadradic hexahedral edgept = np.array([[-1, -1, -1], [ 1, -1, -1], [ 1, 1, -1], @@ -2881,8 +2881,8 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, [ 1, -1, 1], [ 1, 1, 1], [-1, 1, 1]]) - - # midside node indices + + # midside node indices idx = np.array([[0, 1], # 8 [1, 2],# 9 [2, 3],# 10 @@ -2895,14 +2895,14 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, [1, 5],# 17 [2, 6],# 18 [3, 7]], np.int_)# 19 - + midpt = (edgept[idx[:, 0]] + edgept[idx[:, 1]])/2 - + isopts = np.empty((20, 3)) isopts[:8] = edgept isopts[8:] = midpt - - + + # Shape function evaluated at sampling locations u = np.empty((20, 20)) for i in range(20): @@ -2914,8 +2914,8 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, zi**2*eta, zi*eta**2, eta**2*zeta, eta*zeta**2, zeta**2*zi, zeta*zi**2, zi*eta*zeta, zi**2*eta*zeta, zi*eta**2*zeta, zi*eta*zeta**2] - - + + zi, eta, zeta = sympy.symbols('zi eta zeta', real=True) shp = [1, zi, @@ -2937,9 +2937,9 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, eta*zeta*zi**2, eta**2*zeta*zi, eta*zeta**2*zi] - + shape_functions = np.array(np.matrix(shp)*np.matrix(np.linalg.inv(u))).ravel() - + # Take the differential of each shape function evaluated at the # integration point with respect to the natural shape coordinate # system. @@ -2948,8 +2948,8 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, for i in range(20): for j in range(3): nprime[j, i] = sympy.diff(shape_functions[i], variables[j]) - - + + # This is the N_prime functions with the isopts subbed into them pre_j = np.empty((60, 20)) for i in range(20): @@ -2958,7 +2958,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pre_j[i*3 + j, k] = nprime[j, k].subs([(zi, isopts[i, 0]), (eta, isopts[i, 1]), (zeta, isopts[i, 2])]) - + hex_quad_edges = np.empty((20, 3, 8), np.int32) hex_quad_edges[0][0][0] = 1 @@ -2967,61 +2967,61 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[0][1][0] = 8 hex_quad_edges[0][1][1] = 11 hex_quad_edges[0][1][2] = 16 - + hex_quad_edges[1][0][0] = 2 hex_quad_edges[1][0][1] = 0 hex_quad_edges[1][0][2] = 5 hex_quad_edges[1][1][0] = 9 hex_quad_edges[1][1][1] = 8 hex_quad_edges[1][1][2] = 17 - + hex_quad_edges[2][0][0] = 3 hex_quad_edges[2][0][1] = 1 hex_quad_edges[2][0][2] = 6 hex_quad_edges[2][1][0] = 10 hex_quad_edges[2][1][1] = 9 hex_quad_edges[2][1][2] = 18 - + hex_quad_edges[3][0][0] = 0 hex_quad_edges[3][0][1] = 2 hex_quad_edges[3][0][2] = 7 hex_quad_edges[3][1][0] = 11 hex_quad_edges[3][1][1] = 10 hex_quad_edges[3][1][2] = 19 - + hex_quad_edges[4][0][0] = 7 hex_quad_edges[4][0][1] = 5 hex_quad_edges[4][0][2] = 0 hex_quad_edges[4][1][0] = 15 hex_quad_edges[4][1][1] = 12 hex_quad_edges[4][1][2] = 16 - + hex_quad_edges[5][0][0] = 4 hex_quad_edges[5][0][1] = 6 hex_quad_edges[5][0][2] = 1 hex_quad_edges[5][1][0] = 12 hex_quad_edges[5][1][1] = 13 hex_quad_edges[5][1][2] = 17 - + hex_quad_edges[6][0][0] = 5 hex_quad_edges[6][0][1] = 7 hex_quad_edges[6][0][2] = 2 hex_quad_edges[6][1][0] = 13 hex_quad_edges[6][1][1] = 14 hex_quad_edges[6][1][2] = 18 - + hex_quad_edges[7][0][0] = 6 hex_quad_edges[7][0][1] = 4 hex_quad_edges[7][0][2] = 3 hex_quad_edges[7][1][0] = 14 hex_quad_edges[7][1][1] = 15 hex_quad_edges[7][1][2] = 19 - - + + # Node 8 (between edges 0 and 1) hex_quad_edges[8][0][0] = 1 hex_quad_edges[8][0][1] = 0 - + hex_quad_edges[8][1][0] = 9 hex_quad_edges[8][1][1] = 11 hex_quad_edges[8][1][2] = 0 @@ -3030,7 +3030,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[8][1][5] = 3 hex_quad_edges[8][1][6] = 8 hex_quad_edges[8][1][7] = 10 - + hex_quad_edges[8][2][0] = 16 hex_quad_edges[8][2][1] = 17 hex_quad_edges[8][2][2] = 0 @@ -3039,8 +3039,8 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[8][2][5] = 5 hex_quad_edges[8][2][6] = 8 hex_quad_edges[8][2][7] = 12 - - + + # Node 9 (between edges 1 and 2) hex_quad_edges[9][1][0] = 8 hex_quad_edges[9][1][1] = 10 @@ -3050,10 +3050,10 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[9][1][5] = 3 hex_quad_edges[9][1][6] = 9 hex_quad_edges[9][1][7] = 11 - + hex_quad_edges[9][0][0] = 2 hex_quad_edges[9][0][1] = 1 - + hex_quad_edges[9][2][0] = 17 hex_quad_edges[9][2][1] = 18 hex_quad_edges[9][2][2] = 1 @@ -3062,12 +3062,12 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[9][2][5] = 6 hex_quad_edges[9][2][6] = 9 hex_quad_edges[9][2][7] = 13 - - + + # Node 10 (between edges 2 and 3) hex_quad_edges[10][0][0] = 3 hex_quad_edges[10][0][1] = 2 - + hex_quad_edges[10][1][0] = 9 hex_quad_edges[10][1][1] = 11 hex_quad_edges[10][1][2] = 0 @@ -3076,7 +3076,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[10][1][5] = 3 hex_quad_edges[10][1][6] = 10 hex_quad_edges[10][1][7] = 8 - + hex_quad_edges[10][2][0] = 18 hex_quad_edges[10][2][1] = 19 hex_quad_edges[10][2][2] = 2 @@ -3085,13 +3085,13 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[10][2][5] = 7 hex_quad_edges[10][2][6] = 10 hex_quad_edges[10][2][7] = 14 - - - + + + # Node 11 (between edge nodes 0 and 3) hex_quad_edges[11][0][0] = 3 hex_quad_edges[11][0][1] = 0 - + hex_quad_edges[11][1][0] = 16 hex_quad_edges[11][1][1] = 19 hex_quad_edges[11][1][2] = 0 @@ -3100,7 +3100,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[11][1][5] = 7 hex_quad_edges[11][1][6] = 11 hex_quad_edges[11][1][7] = 15 - + hex_quad_edges[11][2][0] = 8 hex_quad_edges[11][2][1] = 10 hex_quad_edges[11][2][2] = 0 @@ -3109,12 +3109,12 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[11][2][5] = 3 hex_quad_edges[11][2][6] = 11 hex_quad_edges[11][2][7] = 9 - - + + # Node 12 (between edges 3 and 0) hex_quad_edges[12][0][0] = 4 hex_quad_edges[12][0][1] = 5 - + hex_quad_edges[12][1][0] = 13 hex_quad_edges[12][1][1] = 15 hex_quad_edges[12][1][2] = 4 @@ -3123,7 +3123,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[12][1][5] = 7 hex_quad_edges[12][1][6] = 12 hex_quad_edges[12][1][7] = 14 - + hex_quad_edges[12][2][0] = 16 hex_quad_edges[12][2][1] = 17 hex_quad_edges[12][2][2] = 0 @@ -3132,8 +3132,8 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[12][2][5] = 5 hex_quad_edges[12][2][6] = 12 hex_quad_edges[12][2][7] = 8 - - + + # Node 13 (between edges 1 and 2) hex_quad_edges[13][1][0] = 12 hex_quad_edges[13][1][1] = 14 @@ -3143,10 +3143,10 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[13][1][5] = 7 hex_quad_edges[13][1][6] = 13 hex_quad_edges[13][1][7] = 15 - + hex_quad_edges[13][0][0] = 5 hex_quad_edges[13][0][1] = 6 - + hex_quad_edges[13][2][0] = 17 hex_quad_edges[13][2][1] = 18 hex_quad_edges[13][2][2] = 1 @@ -3155,12 +3155,12 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[13][2][5] = 6 hex_quad_edges[13][2][6] = 13 hex_quad_edges[13][2][7] = 9 - - + + # Node 14 (between edges 2 and 3) hex_quad_edges[14][0][0] = 6 hex_quad_edges[14][0][1] = 7 - + hex_quad_edges[14][1][0] = 13 hex_quad_edges[14][1][1] = 15 hex_quad_edges[14][1][2] = 4 @@ -3169,7 +3169,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[14][1][5] = 7 hex_quad_edges[14][1][6] = 14 hex_quad_edges[14][1][7] = 12 - + hex_quad_edges[14][2][0] = 18 hex_quad_edges[14][2][1] = 19 hex_quad_edges[14][2][2] = 2 @@ -3178,8 +3178,8 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[14][2][5] = 7 hex_quad_edges[14][2][6] = 14 hex_quad_edges[14][2][7] = 10 - - + + # Node 15 (between edges 1 and 2) hex_quad_edges[15][1][0] = 14 hex_quad_edges[15][1][1] = 12 @@ -3189,10 +3189,10 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[15][1][5] = 7 hex_quad_edges[15][1][6] = 15 hex_quad_edges[15][1][7] = 13 - + hex_quad_edges[15][0][0] = 7 hex_quad_edges[15][0][1] = 4 - + hex_quad_edges[15][2][0] = 16 hex_quad_edges[15][2][1] = 19 hex_quad_edges[15][2][2] = 0 @@ -3201,11 +3201,11 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[15][2][5] = 7 hex_quad_edges[15][2][6] = 15 hex_quad_edges[15][2][7] = 11 - + # Node 16 (between edges 4 and 0) hex_quad_edges[16][0][0] = 0 hex_quad_edges[16][0][1] = 4 - + hex_quad_edges[16][1][0] = 11 hex_quad_edges[16][1][1] = 15 hex_quad_edges[16][1][2] = 0 @@ -3214,7 +3214,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[16][1][5] = 7 hex_quad_edges[16][1][6] = 16 hex_quad_edges[16][1][7] = 19 - + hex_quad_edges[16][2][0] = 8 hex_quad_edges[16][2][1] = 12 hex_quad_edges[16][2][2] = 0 @@ -3223,12 +3223,12 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[16][2][5] = 5 hex_quad_edges[16][2][6] = 16 hex_quad_edges[16][2][7] = 17 - - + + # Node 17 (between edges 4 and 0) hex_quad_edges[17][0][0] = 1 hex_quad_edges[17][0][1] = 5 - + hex_quad_edges[17][1][0] = 8 hex_quad_edges[17][1][1] = 12 hex_quad_edges[17][1][2] = 0 @@ -3237,7 +3237,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[17][1][5] = 5 hex_quad_edges[17][1][6] = 17 hex_quad_edges[17][1][7] = 16 - + hex_quad_edges[17][2][0] = 9 hex_quad_edges[17][2][1] = 13 hex_quad_edges[17][2][2] = 1 @@ -3246,12 +3246,12 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[17][2][5] = 6 hex_quad_edges[17][2][6] = 17 hex_quad_edges[17][2][7] = 18 - - + + # Node 18 (between edges 4 and 0) hex_quad_edges[18][0][0] = 2 hex_quad_edges[18][0][1] = 6 - + hex_quad_edges[18][1][0] = 9 hex_quad_edges[18][1][1] = 13 hex_quad_edges[18][1][2] = 1 @@ -3260,7 +3260,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[18][1][5] = 6 hex_quad_edges[18][1][6] = 18 hex_quad_edges[18][1][7] = 17 - + hex_quad_edges[18][2][0] = 10 hex_quad_edges[18][2][1] = 14 hex_quad_edges[18][2][2] = 2 @@ -3269,12 +3269,12 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[18][2][5] = 7 hex_quad_edges[18][2][6] = 18 hex_quad_edges[18][2][7] = 19 - - + + # Node 19 (between edges 4 and 0) hex_quad_edges[19][0][0] = 3 hex_quad_edges[19][0][1] = 7 - + hex_quad_edges[19][1][0] = 10 hex_quad_edges[19][1][1] = 14 hex_quad_edges[19][1][2] = 2 @@ -3283,7 +3283,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[19][1][5] = 7 hex_quad_edges[19][1][6] = 19 hex_quad_edges[19][1][7] = 18 - + hex_quad_edges[19][2][0] = 11 hex_quad_edges[19][2][1] = 15 hex_quad_edges[19][2][2] = 0 @@ -3292,38 +3292,38 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, hex_quad_edges[19][2][5] = 7 hex_quad_edges[19][2][6] = 19 hex_quad_edges[19][2][7] = 16 - - - + + + for i in range(20): print 'cdef int ind{:d} = cells[c + {:d}]'.format(i, i) - - - + + + for i in range(8): print '# Node {:d}'.format(i) print 'for j in range(3):' print ' e0[j] = 4*pts[ind{:d}, j] - pts[ind{:d}, j] - 3*pts[ind{:d}, j]'.format(hex_quad_edges[i][1][0], hex_quad_edges[i][0][0], i) print ' e1[j] = 4*pts[ind{:d}, j] - pts[ind{:d}, j] - 3*pts[ind{:d}, j]'.format(hex_quad_edges[i][1][1], hex_quad_edges[i][0][1], i) print ' e2[j] = 4*pts[ind{:d}, j] - pts[ind{:d}, j] - 3*pts[ind{:d}, j]'.format(hex_quad_edges[i][1][2], hex_quad_edges[i][0][2], i) - print + print print '# normalize the determinant of the jacobian' - print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' + print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' print 'normjac = triple_product(e1, e2, e0)/tnorm' - print + print print '# Track minimum jacobian' print 'if normjac < jac:' print ' jac = normjac' - print - - + print + + # midside nodes along x moving edges for i in range(8, 20): print '# Node {:d}'.format(i) print 'for j in range(3):' print ' e0[j] = pts[ind{:d}, j] - pts[ind{:d}, j]'.format(hex_quad_edges[i][0][0], hex_quad_edges[i][0][1]) print - + print ' e1[j] = 2*pts[ind{:d}, j] + \\'.format(hex_quad_edges[i][1][0]) print ' 2*pts[ind{:d}, j] - \\'.format(hex_quad_edges[i][1][1]) print ' pts[ind{:d}, j] - \\'.format(hex_quad_edges[i][1][2]) @@ -3341,16 +3341,16 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, print ' pts[ind{:d}, j] - \\'.format(hex_quad_edges[i][2][5]) print ' pts[ind{:d}, j] + \\'.format(hex_quad_edges[i][2][6]) print ' pts[ind{:d}, j]'.format(hex_quad_edges[i][2][7]) - print + print print '# normalize the determinant of the jacobian' - print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' + print 'tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2))' print 'normjac = triple_product(e1, e2, e0)/tnorm' - print + print print '# Track minimum jacobian' print 'if normjac < jac:' print ' jac = normjac' - print - + print + """ cdef double jac = 1.1 @@ -3359,9 +3359,9 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, cdef double [3] e1 cdef double [3] e2 - cdef int i, j + cdef int i, j cdef double normjac, tnorm - + # Store cell indices cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] @@ -3383,123 +3383,123 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, cdef int ind17 = cells[c + 17] cdef int ind18 = cells[c + 18] cdef int ind19 = cells[c + 19] - + # Node 0 for j in range(3): e0[j] = 4*pts[ind8, j] - pts[ind1, j] - 3*pts[ind0, j] e1[j] = 4*pts[ind11, j] - pts[ind3, j] - 3*pts[ind0, j] e2[j] = 4*pts[ind16, j] - pts[ind4, j] - 3*pts[ind0, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 1 for j in range(3): e0[j] = 4*pts[ind9, j] - pts[ind2, j] - 3*pts[ind1, j] e1[j] = 4*pts[ind8, j] - pts[ind0, j] - 3*pts[ind1, j] e2[j] = 4*pts[ind17, j] - pts[ind5, j] - 3*pts[ind1, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 2 for j in range(3): e0[j] = 4*pts[ind10, j] - pts[ind3, j] - 3*pts[ind2, j] e1[j] = 4*pts[ind9, j] - pts[ind1, j] - 3*pts[ind2, j] e2[j] = 4*pts[ind18, j] - pts[ind6, j] - 3*pts[ind2, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 3 for j in range(3): e0[j] = 4*pts[ind11, j] - pts[ind0, j] - 3*pts[ind3, j] e1[j] = 4*pts[ind10, j] - pts[ind2, j] - 3*pts[ind3, j] e2[j] = 4*pts[ind19, j] - pts[ind7, j] - 3*pts[ind3, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 4 for j in range(3): e0[j] = 4*pts[ind15, j] - pts[ind7, j] - 3*pts[ind4, j] e1[j] = 4*pts[ind12, j] - pts[ind5, j] - 3*pts[ind4, j] e2[j] = 4*pts[ind16, j] - pts[ind0, j] - 3*pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 5 for j in range(3): e0[j] = 4*pts[ind12, j] - pts[ind4, j] - 3*pts[ind5, j] e1[j] = 4*pts[ind13, j] - pts[ind6, j] - 3*pts[ind5, j] e2[j] = 4*pts[ind17, j] - pts[ind1, j] - 3*pts[ind5, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 6 for j in range(3): e0[j] = 4*pts[ind13, j] - pts[ind5, j] - 3*pts[ind6, j] e1[j] = 4*pts[ind14, j] - pts[ind7, j] - 3*pts[ind6, j] e2[j] = 4*pts[ind18, j] - pts[ind2, j] - 3*pts[ind6, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 7 for j in range(3): e0[j] = 4*pts[ind14, j] - pts[ind6, j] - 3*pts[ind7, j] e1[j] = 4*pts[ind15, j] - pts[ind4, j] - 3*pts[ind7, j] e2[j] = 4*pts[ind19, j] - pts[ind3, j] - 3*pts[ind7, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 8 for j in range(3): e0[j] = pts[ind1, j] - pts[ind0, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -3508,7 +3508,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind8, j] + \ pts[ind10, j] - + e2[j] = 2*pts[ind16, j] + \ 2*pts[ind17, j] - \ pts[ind0, j] - \ @@ -3517,19 +3517,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind8, j] + \ pts[ind12, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 9 for j in range(3): e0[j] = pts[ind2, j] - pts[ind1, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind10, j] - \ pts[ind0, j] - \ @@ -3538,7 +3538,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind9, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind17, j] + \ 2*pts[ind18, j] - \ pts[ind1, j] - \ @@ -3547,19 +3547,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind9, j] + \ pts[ind13, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 10 for j in range(3): e0[j] = pts[ind3, j] - pts[ind2, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -3568,7 +3568,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind10, j] + \ pts[ind8, j] - + e2[j] = 2*pts[ind18, j] + \ 2*pts[ind19, j] - \ pts[ind2, j] - \ @@ -3577,19 +3577,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind10, j] + \ pts[ind14, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 11 for j in range(3): e0[j] = pts[ind3, j] - pts[ind0, j] - + e1[j] = 2*pts[ind16, j] + \ 2*pts[ind19, j] - \ pts[ind0, j] - \ @@ -3598,7 +3598,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind11, j] + \ pts[ind15, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind10, j] - \ pts[ind0, j] - \ @@ -3607,19 +3607,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind11, j] + \ pts[ind9, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 12 for j in range(3): e0[j] = pts[ind4, j] - pts[ind5, j] - + e1[j] = 2*pts[ind13, j] + \ 2*pts[ind15, j] - \ pts[ind4, j] - \ @@ -3628,7 +3628,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind12, j] + \ pts[ind14, j] - + e2[j] = 2*pts[ind16, j] + \ 2*pts[ind17, j] - \ pts[ind0, j] - \ @@ -3637,19 +3637,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind12, j] + \ pts[ind8, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 13 for j in range(3): e0[j] = pts[ind5, j] - pts[ind6, j] - + e1[j] = 2*pts[ind12, j] + \ 2*pts[ind14, j] - \ pts[ind4, j] - \ @@ -3658,7 +3658,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind13, j] + \ pts[ind15, j] - + e2[j] = 2*pts[ind17, j] + \ 2*pts[ind18, j] - \ pts[ind1, j] - \ @@ -3667,19 +3667,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind13, j] + \ pts[ind9, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 14 for j in range(3): e0[j] = pts[ind6, j] - pts[ind7, j] - + e1[j] = 2*pts[ind13, j] + \ 2*pts[ind15, j] - \ pts[ind4, j] - \ @@ -3688,7 +3688,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind14, j] + \ pts[ind12, j] - + e2[j] = 2*pts[ind18, j] + \ 2*pts[ind19, j] - \ pts[ind2, j] - \ @@ -3697,19 +3697,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind14, j] + \ pts[ind10, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 15 for j in range(3): e0[j] = pts[ind7, j] - pts[ind4, j] - + e1[j] = 2*pts[ind14, j] + \ 2*pts[ind12, j] - \ pts[ind4, j] - \ @@ -3718,7 +3718,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind15, j] + \ pts[ind13, j] - + e2[j] = 2*pts[ind16, j] + \ 2*pts[ind19, j] - \ pts[ind0, j] - \ @@ -3727,19 +3727,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind15, j] + \ pts[ind11, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 16 for j in range(3): e0[j] = pts[ind0, j] - pts[ind4, j] - + e1[j] = 2*pts[ind11, j] + \ 2*pts[ind15, j] - \ pts[ind0, j] - \ @@ -3748,7 +3748,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind16, j] + \ pts[ind19, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind12, j] - \ pts[ind0, j] - \ @@ -3757,19 +3757,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind16, j] + \ pts[ind17, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 17 for j in range(3): e0[j] = pts[ind1, j] - pts[ind5, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind12, j] - \ pts[ind0, j] - \ @@ -3778,7 +3778,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind17, j] + \ pts[ind16, j] - + e2[j] = 2*pts[ind9, j] + \ 2*pts[ind13, j] - \ pts[ind1, j] - \ @@ -3787,19 +3787,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind17, j] + \ pts[ind18, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 18 for j in range(3): e0[j] = pts[ind2, j] - pts[ind6, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind13, j] - \ pts[ind1, j] - \ @@ -3808,7 +3808,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind18, j] + \ pts[ind17, j] - + e2[j] = 2*pts[ind10, j] + \ 2*pts[ind14, j] - \ pts[ind2, j] - \ @@ -3817,19 +3817,19 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind18, j] + \ pts[ind19, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 19 for j in range(3): e0[j] = pts[ind3, j] - pts[ind7, j] - + e1[j] = 2*pts[ind10, j] + \ 2*pts[ind14, j] - \ pts[ind2, j] - \ @@ -3838,7 +3838,7 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind19, j] + \ pts[ind18, j] - + e2[j] = 2*pts[ind11, j] + \ 2*pts[ind15, j] - \ pts[ind0, j] - \ @@ -3847,15 +3847,15 @@ cdef inline double hex_quad_qual(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind19, j] + \ pts[ind16, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm(e0)*vector_norm(e1)*vector_norm(e2)) normjac = triple_product(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + return jac @@ -3865,9 +3865,9 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, cdef float [3] e0 cdef float [3] e1 cdef float [3] e2 - cdef int i, j + cdef int i, j cdef float normjac, tnorm - + # Store cell indices cdef int ind0 = cells[c + 0] cdef int ind1 = cells[c + 1] @@ -3889,123 +3889,123 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, cdef int ind17 = cells[c + 17] cdef int ind18 = cells[c + 18] cdef int ind19 = cells[c + 19] - + # Node 0 for j in range(3): e0[j] = 4*pts[ind8, j] - pts[ind1, j] - 3*pts[ind0, j] e1[j] = 4*pts[ind11, j] - pts[ind3, j] - 3*pts[ind0, j] e2[j] = 4*pts[ind16, j] - pts[ind4, j] - 3*pts[ind0, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 1 for j in range(3): e0[j] = 4*pts[ind9, j] - pts[ind2, j] - 3*pts[ind1, j] e1[j] = 4*pts[ind8, j] - pts[ind0, j] - 3*pts[ind1, j] e2[j] = 4*pts[ind17, j] - pts[ind5, j] - 3*pts[ind1, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 2 for j in range(3): e0[j] = 4*pts[ind10, j] - pts[ind3, j] - 3*pts[ind2, j] e1[j] = 4*pts[ind9, j] - pts[ind1, j] - 3*pts[ind2, j] e2[j] = 4*pts[ind18, j] - pts[ind6, j] - 3*pts[ind2, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 3 for j in range(3): e0[j] = 4*pts[ind11, j] - pts[ind0, j] - 3*pts[ind3, j] e1[j] = 4*pts[ind10, j] - pts[ind2, j] - 3*pts[ind3, j] e2[j] = 4*pts[ind19, j] - pts[ind7, j] - 3*pts[ind3, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 4 for j in range(3): e0[j] = 4*pts[ind15, j] - pts[ind7, j] - 3*pts[ind4, j] e1[j] = 4*pts[ind12, j] - pts[ind5, j] - 3*pts[ind4, j] e2[j] = 4*pts[ind16, j] - pts[ind0, j] - 3*pts[ind4, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 5 for j in range(3): e0[j] = 4*pts[ind12, j] - pts[ind4, j] - 3*pts[ind5, j] e1[j] = 4*pts[ind13, j] - pts[ind6, j] - 3*pts[ind5, j] e2[j] = 4*pts[ind17, j] - pts[ind1, j] - 3*pts[ind5, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 6 for j in range(3): e0[j] = 4*pts[ind13, j] - pts[ind5, j] - 3*pts[ind6, j] e1[j] = 4*pts[ind14, j] - pts[ind7, j] - 3*pts[ind6, j] e2[j] = 4*pts[ind18, j] - pts[ind2, j] - 3*pts[ind6, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 7 for j in range(3): e0[j] = 4*pts[ind14, j] - pts[ind6, j] - 3*pts[ind7, j] e1[j] = 4*pts[ind15, j] - pts[ind4, j] - 3*pts[ind7, j] e2[j] = 4*pts[ind19, j] - pts[ind3, j] - 3*pts[ind7, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 8 for j in range(3): e0[j] = pts[ind1, j] - pts[ind0, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -4014,7 +4014,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind8, j] + \ pts[ind10, j] - + e2[j] = 2*pts[ind16, j] + \ 2*pts[ind17, j] - \ pts[ind0, j] - \ @@ -4023,19 +4023,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind8, j] + \ pts[ind12, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 9 for j in range(3): e0[j] = pts[ind2, j] - pts[ind1, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind10, j] - \ pts[ind0, j] - \ @@ -4044,7 +4044,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind9, j] + \ pts[ind11, j] - + e2[j] = 2*pts[ind17, j] + \ 2*pts[ind18, j] - \ pts[ind1, j] - \ @@ -4053,19 +4053,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind9, j] + \ pts[ind13, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 10 for j in range(3): e0[j] = pts[ind3, j] - pts[ind2, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind11, j] - \ pts[ind0, j] - \ @@ -4074,7 +4074,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind10, j] + \ pts[ind8, j] - + e2[j] = 2*pts[ind18, j] + \ 2*pts[ind19, j] - \ pts[ind2, j] - \ @@ -4083,19 +4083,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind10, j] + \ pts[ind14, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 11 for j in range(3): e0[j] = pts[ind3, j] - pts[ind0, j] - + e1[j] = 2*pts[ind16, j] + \ 2*pts[ind19, j] - \ pts[ind0, j] - \ @@ -4104,7 +4104,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind11, j] + \ pts[ind15, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind10, j] - \ pts[ind0, j] - \ @@ -4113,19 +4113,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind3, j] - \ pts[ind11, j] + \ pts[ind9, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 12 for j in range(3): e0[j] = pts[ind4, j] - pts[ind5, j] - + e1[j] = 2*pts[ind13, j] + \ 2*pts[ind15, j] - \ pts[ind4, j] - \ @@ -4134,7 +4134,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind12, j] + \ pts[ind14, j] - + e2[j] = 2*pts[ind16, j] + \ 2*pts[ind17, j] - \ pts[ind0, j] - \ @@ -4143,19 +4143,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind12, j] + \ pts[ind8, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 13 for j in range(3): e0[j] = pts[ind5, j] - pts[ind6, j] - + e1[j] = 2*pts[ind12, j] + \ 2*pts[ind14, j] - \ pts[ind4, j] - \ @@ -4164,7 +4164,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind13, j] + \ pts[ind15, j] - + e2[j] = 2*pts[ind17, j] + \ 2*pts[ind18, j] - \ pts[ind1, j] - \ @@ -4173,19 +4173,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind13, j] + \ pts[ind9, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 14 for j in range(3): e0[j] = pts[ind6, j] - pts[ind7, j] - + e1[j] = 2*pts[ind13, j] + \ 2*pts[ind15, j] - \ pts[ind4, j] - \ @@ -4194,7 +4194,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind14, j] + \ pts[ind12, j] - + e2[j] = 2*pts[ind18, j] + \ 2*pts[ind19, j] - \ pts[ind2, j] - \ @@ -4203,19 +4203,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind14, j] + \ pts[ind10, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 15 for j in range(3): e0[j] = pts[ind7, j] - pts[ind4, j] - + e1[j] = 2*pts[ind14, j] + \ 2*pts[ind12, j] - \ pts[ind4, j] - \ @@ -4224,7 +4224,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind15, j] + \ pts[ind13, j] - + e2[j] = 2*pts[ind16, j] + \ 2*pts[ind19, j] - \ pts[ind0, j] - \ @@ -4233,19 +4233,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind15, j] + \ pts[ind11, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 16 for j in range(3): e0[j] = pts[ind0, j] - pts[ind4, j] - + e1[j] = 2*pts[ind11, j] + \ 2*pts[ind15, j] - \ pts[ind0, j] - \ @@ -4254,7 +4254,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind16, j] + \ pts[ind19, j] - + e2[j] = 2*pts[ind8, j] + \ 2*pts[ind12, j] - \ pts[ind0, j] - \ @@ -4263,19 +4263,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind16, j] + \ pts[ind17, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 17 for j in range(3): e0[j] = pts[ind1, j] - pts[ind5, j] - + e1[j] = 2*pts[ind8, j] + \ 2*pts[ind12, j] - \ pts[ind0, j] - \ @@ -4284,7 +4284,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind5, j] - \ pts[ind17, j] + \ pts[ind16, j] - + e2[j] = 2*pts[ind9, j] + \ 2*pts[ind13, j] - \ pts[ind1, j] - \ @@ -4293,19 +4293,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind17, j] + \ pts[ind18, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 18 for j in range(3): e0[j] = pts[ind2, j] - pts[ind6, j] - + e1[j] = 2*pts[ind9, j] + \ 2*pts[ind13, j] - \ pts[ind1, j] - \ @@ -4314,7 +4314,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind6, j] - \ pts[ind18, j] + \ pts[ind17, j] - + e2[j] = 2*pts[ind10, j] + \ 2*pts[ind14, j] - \ pts[ind2, j] - \ @@ -4323,19 +4323,19 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind18, j] + \ pts[ind19, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + # Node 19 for j in range(3): e0[j] = pts[ind3, j] - pts[ind7, j] - + e1[j] = 2*pts[ind10, j] + \ 2*pts[ind14, j] - \ pts[ind2, j] - \ @@ -4344,7 +4344,7 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind19, j] + \ pts[ind18, j] - + e2[j] = 2*pts[ind11, j] + \ 2*pts[ind15, j] - \ pts[ind0, j] - \ @@ -4353,15 +4353,15 @@ cdef inline float hex_quad_qual_float(int64_t [::1] cells, int c, pts[ind7, j] - \ pts[ind19, j] + \ pts[ind16, j] - + # normalize the determinant of the jacobian tnorm = (vector_norm_float(e0)*vector_norm_float(e1)*vector_norm_float(e2)) normjac = triple_product_float(e1, e2, e0)/tnorm - + # Track minimum jacobian if normjac < jac: jac = normjac - + return jac @@ -4373,12 +4373,12 @@ cdef inline double volume_weg(double [8][3] cell_points) nogil: vector_sub(cell_points[1], cell_points[2], e1) vector_sub(cell_points[2], cell_points[4], e2) cdef double volume = triple_product(e0, e1, e2) - + vector_sub(cell_points[0], cell_points[4], e0) vector_sub(cell_points[4], cell_points[2], e1) vector_sub(cell_points[2], cell_points[5], e2) volume += triple_product(e0, e1, e2) - + vector_sub(cell_points[4], cell_points[5], e1) vector_sub(cell_points[5], cell_points[3], e2) volume += triple_product(e0, e1, e2) @@ -4393,12 +4393,12 @@ cdef inline float volume_weg_float(float [8][3] cell_points) nogil: vector_sub_float(cell_points[1], cell_points[2], e1) vector_sub_float(cell_points[2], cell_points[4], e2) cdef float volume = triple_product_float(e0, e1, e2) - + vector_sub_float(cell_points[0], cell_points[4], e0) vector_sub_float(cell_points[4], cell_points[2], e1) vector_sub_float(cell_points[2], cell_points[5], e2) volume += triple_product_float(e0, e1, e2) - + vector_sub_float(cell_points[4], cell_points[5], e1) vector_sub_float(cell_points[5], cell_points[3], e2) volume += triple_product_float(e0, e1, e2) @@ -4479,10 +4479,10 @@ cdef inline float volume_hex_float(float [8][3] points) nogil: def cell_quality(int64_t [::1] cells, int64_t [::1] offset, uint8 [::1] celltypes, double [:, ::1] pts, int as_linear=0): - """ + """ Returns the minimum scaled jacobian for each cell given a cell array from a vtk unstructured grid. Accounts for midside nodes. - + Parameters ---------- cells : int64_t [::1] @@ -4502,7 +4502,7 @@ def cell_quality(int64_t [::1] cells, int64_t [::1] offset, ------- quality : np.ndarray Minimum scaled jacobian for each valid cell. - + Notes ----- Checks the quality of the following cell types: @@ -4569,7 +4569,7 @@ def cell_quality(int64_t [::1] cells, int64_t [::1] offset, def cell_quality_float(int64_t [::1] cells, int64_t [::1] offset, uint8 [::1] celltypes, float [:, ::1] pts, int as_linear=0): - """ + """ Returns the minimum scaled jacobian for each cell given a cell array from a vtk unstructured grid. Accounts for midside nodes. @@ -4592,7 +4592,7 @@ def cell_quality_float(int64_t [::1] cells, int64_t [::1] offset, ------- quality : np.ndarray Minimum scaled jacobian for each valid cell. - + Notes ----- Checks the quality of the following cell types: diff --git a/ansys/mapdl/reader/cython/_parsefull.pyx b/ansys/mapdl/reader/cython/_parsefull.pyx index f961d0ed..b44ac585 100644 --- a/ansys/mapdl/reader/cython/_parsefull.pyx +++ b/ansys/mapdl/reader/cython/_parsefull.pyx @@ -12,7 +12,7 @@ import numpy as np # Definitions from c header cdef extern from "parsefull.h": int return_fheader(char*, int*) - + void read_full(int*, int*, int*, int*, int*, double*, int*, int*, double*, int*, char*, int*, int); @@ -22,25 +22,25 @@ def return_header(filename): # Convert python string to char array cdef bytes py_bytes = filename.encode() cdef char* c_filename = py_bytes - + # Read header cdef int [::1] fheader = np.empty(101, np.int32) return_fheader(c_filename, &fheader[0]) - + return np.asarray(fheader) - - + + def load_km(filename, is_sorted): - """Reads an ANSYS full file and returns indices to construct symmetric, real, + """Reads an ANSYS full file and returns indices to construct symmetric, real, and sparse mass and stiffness matrices """ # convert to int cdef int sort = is_sorted - + # Convert python string to char array cdef bytes py_bytes = filename.encode() cdef char* c_filename = py_bytes - + # Read header cdef int [::1] fheader = np.empty(101, np.int32) cdef int rst = return_fheader(c_filename, &fheader[0]) @@ -89,7 +89,7 @@ def load_km(filename, is_sorted): nfree = numdat[0] kentry = numdat[1] mentry = numdat[2] - + # Resort degree of freedom references if sorted cdef int i, sind cdef int[::1] nref_sort = np.empty(nfree, np.int32) @@ -103,17 +103,17 @@ def load_km(filename, is_sorted): # Return sorted arrays to python return (np.asarray(nref_sort), np.asarray(dref_sort), - np.asarray(krows)[:kentry], + np.asarray(krows)[:kentry], np.asarray(kcols)[:kentry], np.asarray(kdata)[:kentry], np.asarray(mrows)[:mentry], np.asarray(mcols)[:mentry], - np.asarray(mdata)[:mentry]) + np.asarray(mdata)[:mentry]) else: # Return unsorted arrays to python return (np.asarray(nref)[:nfree], np.asarray(dref)[:nfree], - np.asarray(krows)[:kentry], + np.asarray(krows)[:kentry], np.asarray(kcols)[:kentry], np.asarray(kdata)[:kentry], np.asarray(mrows)[:mentry], diff --git a/ansys/mapdl/reader/cython/_reader.pyx b/ansys/mapdl/reader/cython/_reader.pyx index 4c7f3cb3..8ec305e0 100644 --- a/ansys/mapdl/reader/cython/_reader.pyx +++ b/ansys/mapdl/reader/cython/_reader.pyx @@ -60,9 +60,9 @@ cdef extern from 'vtk_support.h': cdef int myfgets(char *outstr, char *instr, int64_t *n, int64_t fsize): """Copies a single line from instr to outstr starting from position n """ - + cdef int k = n[0] - + # Search line at a maximum of 10000 characters cdef int64_t i, c c = n[0] @@ -70,7 +70,7 @@ cdef int myfgets(char *outstr, char *instr, int64_t *n, int64_t fsize): # check if end of file if c > fsize: return 1 - + # Add null character if at end of line if instr[c] == '\r': n[0] += i + 2 @@ -80,11 +80,11 @@ cdef int myfgets(char *outstr, char *instr, int64_t *n, int64_t fsize): n[0] += i + 1 outstr[i] = '\0' return 0 - + # Otherwise, store data to output string outstr[i] = instr[c] c += 1 - + # Line exceeds 1000 char (unlikely with ANSYS CDB formatting) return 1 @@ -168,7 +168,7 @@ def read(filename, read_parameters=False, debug=False, read_eblock=True): # File counter cdef int tmpval, start_pos cdef int64_t n = 0 - + # Define variables cdef size_t l = 0 cdef ssize_t read @@ -179,7 +179,7 @@ def read(filename, read_parameters=False, debug=False, read_eblock=True): # Size temp char array cdef char line[1000] cdef char tempstr[100] - + # Get element types elem_type = [] rnum = [] @@ -274,7 +274,7 @@ def read(filename, read_parameters=False, debug=False, read_eblock=True): if debug: print('reading RLBLOCK') - + # The RLBLOCK command defines a real constant # set. The real constant sets follow each set, # starting with @@ -326,10 +326,10 @@ def read(filename, read_parameters=False, debug=False, read_eblock=True): for i in range(6): rcon.append(float(line[16 + 16*i:32 + 16*i])) ncon -= 1 - + # advance line if myfgets(line, raw, &n, fsize): raise Exception(badstr) - + # read next line while True: if ncon > 7: @@ -338,22 +338,22 @@ def read(filename, read_parameters=False, debug=False, read_eblock=True): ncon -= 1 # advance if myfgets(line, raw, &n, fsize): raise Exception(badstr) - + else: for i in range(ncon): - try: - rcon.append(float(line[16*i:16 + 16*i])) + try: + rcon.append(float(line[16*i:16 + 16*i])) # account for empty 0 values except: rcon.append(0.0) - + break - + # If only one in constant data else: for i in range(ncon): - rcon.append(float(line[16 + 16*i:32 + 16*i])) - + rcon.append(float(line[16 + 16*i:32 + 16*i])) + rdat.append(rcon) elif 'N' == line[0] or 'n' == line[0]: @@ -373,7 +373,7 @@ def read(filename, read_parameters=False, debug=False, read_eblock=True): nodes_read = True # Get size of NBLOCK nnodes = int(line[line.rfind(b',') + 1:]) - # this value may be wrong... + # this value may be wrong... # Get format of NBLOCK if myfgets(line, raw, &n, fsize): @@ -549,7 +549,7 @@ def node_block_format(string): nexp = 2 # default when missing nfields = 6 f_size = 21 - c = 0 + c = 0 for field in fields: if 'i' in field: items = field.split('i') @@ -582,7 +582,7 @@ def component_interperter(component): f_new.append(component[i]) else: # otherwise, append list f_new.append(range(abs(component[i - 1]) + 1, abs(component[i]) + 1)) - + return np.hstack(f_new).astype(ctypes.c_int) diff --git a/ansys/mapdl/reader/cython/archive.c b/ansys/mapdl/reader/cython/archive.c index cb2f76a0..ba4bb24e 100644 --- a/ansys/mapdl/reader/cython/archive.c +++ b/ansys/mapdl/reader/cython/archive.c @@ -1,10 +1,10 @@ +#include #include #include -#include // necessary for ubuntu build on azure #ifdef __linux__ - #include +#include #endif // VTK cell types @@ -23,7 +23,7 @@ // Write node IDs, node coordinates, and angles to file as a NBLOCK int write_nblock(FILE *file, const int n_nodes, const int max_node_id, const int *node_id, const double *nodes, const double *angles, - int has_angles){ + int has_angles) { // Header // Tell ANSYS to start reading the node block with 6 fields, // associated with a solid, the maximum node number and the number @@ -34,18 +34,16 @@ int write_nblock(FILE *file, const int n_nodes, const int max_node_id, int i; if (has_angles) { - for (i=0; i +#include +#include +#include #include -#include #include -#include -#include -#include +#include // necessary for ubuntu build on azure #ifdef __linux__ - #include +#include #endif using namespace std; -#define MEM_ZERO(where,size) memset((where),'\0',(size)) -#define MEM_COPY(from,to,size) memcpy((to),(from),(size)) -#define MEMCOPY(from,to,n_items,type) MEM_COPY((char *)(from),(char *)(to),(unsigned)(n_items)*sizeof(type)) -#define IS_ON(e,p) ((e) & (1u << (p))) - - -static int NbBitsOn( int iVal) -{ - static int nbbitsperchar[] = - {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4, - 3,4,4,5,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5, - 3,4,4,5,4,5,5,6,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4, - 3,4,4,5,3,4,4,5,4,5,5,6,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6, - 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,1,2,2,3,2,3,3,4,2,3,3,4, - 3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,2,3,3,4,3,4,4,5, - 3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,2,3,3,4, - 3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7, - 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8}; - - unsigned char *cval = (unsigned char *)(&iVal); - - return( nbbitsperchar[cval[0]] + nbbitsperchar[cval[1]] + - nbbitsperchar[cval[2]] + nbbitsperchar[cval[3]]); +#define MEM_ZERO(where, size) memset((where), '\0', (size)) +#define MEM_COPY(from, to, size) memcpy((to), (from), (size)) +#define MEMCOPY(from, to, n_items, type) \ + MEM_COPY((char *)(from), (char *)(to), (unsigned)(n_items) * sizeof(type)) +#define IS_ON(e, p) ((e) & (1u << (p))) + +static int NbBitsOn(int iVal) { + static int nbbitsperchar[] = { + 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4, 1, 2, 2, 3, 2, 3, 3, 4, + 2, 3, 3, 4, 3, 4, 4, 5, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, + 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 1, 2, 2, 3, 2, 3, 3, 4, + 2, 3, 3, 4, 3, 4, 4, 5, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, + 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, + 4, 5, 5, 6, 5, 6, 6, 7, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5, + 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 2, 3, 3, 4, 3, 4, 4, 5, + 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, + 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6, 3, 4, 4, 5, 4, 5, 5, 6, + 4, 5, 5, 6, 5, 6, 6, 7, 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7, + 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8}; + + unsigned char *cval = (unsigned char *)(&iVal); + + return (nbbitsperchar[cval[0]] + nbbitsperchar[cval[1]] + + nbbitsperchar[cval[2]] + nbbitsperchar[cval[3]]); } - // Read in record and determine size // bsparse_flag true when record uses binary compression // type_flag true when using integers // prec_flag true when using single precision (short for int) -int read_header(ifstream* binFile, int* bsparse_flag, int* wsparse_flag, - int* zlib_flag, int* prec_flag, int* type_flag){ +int read_header(ifstream *binFile, int *bsparse_flag, int *wsparse_flag, + int *zlib_flag, int *prec_flag, int *type_flag) { char *raw = new char[8]; // read the first 8 bytes, includes total buffer size and flags binFile->read(raw, 8); - int bufsize = *(int*)&raw[0]; + int bufsize = *(int *)&raw[0]; // bsparse flag *bsparse_flag = (raw[7] >> 3) & 1; @@ -75,14 +75,14 @@ int read_header(ifstream* binFile, int* bsparse_flag, int* wsparse_flag, // bsparse_flag true when record uses binary compression // type_flag true when using integers // prec_flag true when using single precision (short for int) -int read_header(fstream* binFile, int* bsparse_flag, int* wsparse_flag, - int* zlib_flag, int* prec_flag, int* type_flag){ +int read_header(fstream *binFile, int *bsparse_flag, int *wsparse_flag, + int *zlib_flag, int *prec_flag, int *type_flag) { char *raw = new char[8]; // read the first 8 bytes, includes total buffer size and flags binFile->read(raw, 8); - int bufsize = *(int*)&raw[0]; + int bufsize = *(int *)&raw[0]; // bsparse flag *bsparse_flag = (raw[7] >> 3) & 1; @@ -103,544 +103,530 @@ int read_header(fstream* binFile, int* bsparse_flag, int* wsparse_flag, return bufsize; } - - -template -char* ReadBsparseRecord(T *buffer, int *size){ - int *raw = (int*)buffer; +template char *ReadBsparseRecord(T *buffer, int *size) { + int *raw = (int *)buffer; *size = *raw++; int bitcod = *raw++; - T *tbuf = (T*)raw; + T *tbuf = (T *)raw; T *vec = new T[*size]; int iloc = -1; - while (++iloc < *size){ - if (IS_ON(bitcod, iloc)){ // store value - vec[iloc] = *tbuf++; - } else{ // set value to zero + while (++iloc < *size) { + if (IS_ON(bitcod, iloc)) { // store value + vec[iloc] = *tbuf++; + } else { // set value to zero vec[iloc] = 0; } } - return (char*)vec; + return (char *)vec; } - // read binary sparse record and store in a vector -template -void ReadBsparseRecordToVec(int *raw, int *size, T*vec){ +template void ReadBsparseRecordToVec(int *raw, int *size, T *vec) { *size = *raw++; int bitcod = *raw++; - T *tbuf = (T*)raw; + T *tbuf = (T *)raw; int iloc = -1; - while (++iloc < *size){ - if (IS_ON(bitcod, iloc)){ // store value - vec[iloc] = *tbuf++; - } else{ // set value to zero + while (++iloc < *size) { + if (IS_ON(bitcod, iloc)) { // store value + vec[iloc] = *tbuf++; + } else { // set value to zero vec[iloc] = 0; } } } - - -char* ReadShortBsparseRecord(int *raw, int *size){ +char *ReadShortBsparseRecord(int *raw, int *size) { *size = *raw++; int bitcod = *raw++; short *vec = new short[*size](); - short *tbuf = (short *)raw; + short *tbuf = (short *)raw; int iloc = -1; int nb = NbBitsOn(bitcod); - if (nb%2) nb++; - while ( ++iloc < *size) - { - if ( IS_ON( bitcod, iloc)){ - vec[iloc] = *tbuf++; - } else{ - vec[iloc] = 0; - } + if (nb % 2) + nb++; + while (++iloc < *size) { + if (IS_ON(bitcod, iloc)) { + vec[iloc] = *tbuf++; + } else { + vec[iloc] = 0; } + } - return (char*)vec; + return (char *)vec; } - -void ReadShortBsparseRecordToVec(int *raw, int *size, short *vec) -{ +void ReadShortBsparseRecordToVec(int *raw, int *size, short *vec) { *size = *raw++; int bitcod = *raw++; - short *tbuf = (short *)raw; - int iloc = -1; + short *tbuf = (short *)raw; + int iloc = -1; int nb = NbBitsOn(bitcod); - if (nb%2) nb++; - while ( ++iloc < *size) - { - if ( IS_ON( bitcod, iloc)){ - vec[iloc] = *tbuf++; - } else { - vec[iloc] = 0; - } + if (nb % 2) + nb++; + while (++iloc < *size) { + if (IS_ON(bitcod, iloc)) { + vec[iloc] = *tbuf++; + } else { + vec[iloc] = 0; } - + } } +template char *ReadWindowedSparseBuffer(T *buffer, int *size) { -template -char* ReadWindowedSparseBuffer(T *buffer, int *size){ + int iShift = sizeof(T) / sizeof(int); - int iShift = sizeof(T)/sizeof(int); - - int *raw = (int*)buffer; + int *raw = (int *)buffer; *size = *raw++; int NWin = *raw++; T *vec = new T[*size]; T *adr = vec; - MEM_ZERO( adr, *size*sizeof(T)); + MEM_ZERO(adr, *size * sizeof(T)); int iLoc, iLen; - if ( NWin > 0) - do { + if (NWin > 0) + do { - /* ===== We read the location of the new Window */ - iLoc = *raw++; /* ===== iLoc = Where start the next window */ - if ( iLoc > 0) /* ===== One isolated NonZero Value - no need to store a Window Len */ + /* ===== We read the location of the new Window */ + iLoc = *raw++; /* ===== iLoc = Where start the next window */ + if (iLoc > 0) /* ===== One isolated NonZero Value - no need to store a + Window Len */ { - vec[iLoc] = ((T *)raw)[0]; - raw += iShift; - } - else /* ===== New Window of size iLen */ + vec[iLoc] = ((T *)raw)[0]; + raw += iShift; + } else /* ===== New Window of size iLen */ { - T *adr = vec + (-iLoc); /* Start of the Windows in the output vector */ - iLen = *raw++; /* Length of the Window */ - - if (iLen > 0) /* ===== Non Constant Values */ - { - MEMCOPY( (T *)raw, adr, iLen, T); - - raw += iShift*iLen; - } - else /* ===== Constant Value : only one value is stored */ - { - iLen = - iLen; - T ValCst = *((T *)(raw)); raw += iShift; - do (*(adr++) = ValCst); while (--iLen > 0); - } + T *adr = vec + (-iLoc); /* Start of the Windows in the output vector */ + iLen = *raw++; /* Length of the Window */ + + if (iLen > 0) /* ===== Non Constant Values */ + { + MEMCOPY((T *)raw, adr, iLen, T); + + raw += iShift * iLen; + } else /* ===== Constant Value : only one value is stored */ + { + iLen = -iLen; + T ValCst = *((T *)(raw)); + raw += iShift; + do + (*(adr++) = ValCst); + while (--iLen > 0); + } } - } while ( --NWin > 0); + } while (--NWin > 0); - return (char*)vec; + return (char *)vec; } - -char* ReadWindowedSparseBufferDouble(int *raw, int *size, double *vec){ - int iShift = sizeof(double)/sizeof(int); +char *ReadWindowedSparseBufferDouble(int *raw, int *size, double *vec) { + int iShift = sizeof(double) / sizeof(int); *size = *raw++; int NWin = *raw++; double *adr = vec; - MEM_ZERO( adr, *size*sizeof(double)); + MEM_ZERO(adr, *size * sizeof(double)); int iLoc, iLen; - if ( NWin > 0) - do { + if (NWin > 0) + do { - /* ===== We read the location of the new Window */ - iLoc = *raw++; /* ===== iLoc = Where start the next window */ + /* ===== We read the location of the new Window */ + iLoc = *raw++; /* ===== iLoc = Where start the next window */ - if ( iLoc > 0) /* ===== One isolated NonZero Value - no need to store a Window Len */ + if (iLoc > 0) /* ===== One isolated NonZero Value - no need to store a + Window Len */ { - vec[iLoc] = ((double *)raw)[0]; - raw += iShift; - } - else /* ===== New Window of size iLen */ + vec[iLoc] = ((double *)raw)[0]; + raw += iShift; + } else /* ===== New Window of size iLen */ { - double *adr = vec + (-iLoc); /* Start of the Windows in the output vector */ - iLen = *raw++; /* Length of the Window */ - - if (iLen > 0) /* ===== Non Constant Values */ - { - MEMCOPY( (double *)raw, adr, iLen, double); - - raw += iShift*iLen; - } - else /* ===== Constant Value : only one value is stored */ - { - iLen = - iLen; - double ValCst = *((double *)(raw)); raw += iShift; - do (*(adr++) = ValCst); while (--iLen > 0); - } + double *adr = + vec + (-iLoc); /* Start of the Windows in the output vector */ + iLen = *raw++; /* Length of the Window */ + + if (iLen > 0) /* ===== Non Constant Values */ + { + MEMCOPY((double *)raw, adr, iLen, double); + + raw += iShift * iLen; + } else /* ===== Constant Value : only one value is stored */ + { + iLen = -iLen; + double ValCst = *((double *)(raw)); + raw += iShift; + do + (*(adr++) = ValCst); + while (--iLen > 0); + } } - } while ( --NWin > 0); + } while (--NWin > 0); - return (char*)vec; + return (char *)vec; } - -char* ReadWindowedSparseBufferFloat(int *raw, int *size, float *vec){ - int iShift = sizeof(float)/sizeof(int); +char *ReadWindowedSparseBufferFloat(int *raw, int *size, float *vec) { + int iShift = sizeof(float) / sizeof(int); *size = *raw++; int NWin = *raw++; float *adr = vec; - MEM_ZERO( adr, *size*sizeof(float)); + MEM_ZERO(adr, *size * sizeof(float)); int iLoc, iLen; - if ( NWin > 0) - do { + if (NWin > 0) + do { - /* ===== We read the location of the new Window */ - iLoc = *raw++; /* ===== iLoc = Where start the next window */ + /* ===== We read the location of the new Window */ + iLoc = *raw++; /* ===== iLoc = Where start the next window */ - if ( iLoc > 0) /* ===== One isolated NonZero Value - no need to store a Window Len */ + if (iLoc > 0) /* ===== One isolated NonZero Value - no need to store a + Window Len */ { - vec[iLoc] = ((float *)raw)[0]; - raw += iShift; - } - else /* ===== New Window of size iLen */ + vec[iLoc] = ((float *)raw)[0]; + raw += iShift; + } else /* ===== New Window of size iLen */ { - float *adr = vec + (-iLoc); /* Start of the Windows in the output vector */ - iLen = *raw++; /* Length of the Window */ - - if (iLen > 0) /* ===== Non Constant Values */ - { - MEMCOPY( (float *)raw, adr, iLen, float); - - raw += iShift*iLen; - } - else /* ===== Constant Value : only one value is stored */ - { - iLen = - iLen; - float ValCst = *((float *)(raw)); raw += iShift; - do (*(adr++) = ValCst); while (--iLen > 0); - } + float *adr = + vec + (-iLoc); /* Start of the Windows in the output vector */ + iLen = *raw++; /* Length of the Window */ + + if (iLen > 0) /* ===== Non Constant Values */ + { + MEMCOPY((float *)raw, adr, iLen, float); + + raw += iShift * iLen; + } else /* ===== Constant Value : only one value is stored */ + { + iLen = -iLen; + float ValCst = *((float *)(raw)); + raw += iShift; + do + (*(adr++) = ValCst); + while (--iLen > 0); + } } - } while ( --NWin > 0); + } while (--NWin > 0); - return (char*)vec; + return (char *)vec; } +char *ReadWindowedSparseBufferInt(int *raw, int *size, int *vec) { -char* ReadWindowedSparseBufferInt(int *raw, int *size, int *vec){ - - int iShift = sizeof(int)/sizeof(int); + int iShift = sizeof(int) / sizeof(int); *size = *raw++; int NWin = *raw++; int *adr = vec; - MEM_ZERO( adr, *size*sizeof(int)); + MEM_ZERO(adr, *size * sizeof(int)); int iLoc, iLen; - if ( NWin > 0) - do { + if (NWin > 0) + do { - /* ===== We read the location of the new Window */ - iLoc = *raw++; /* ===== iLoc = Where start the next window */ + /* ===== We read the location of the new Window */ + iLoc = *raw++; /* ===== iLoc = Where start the next window */ - if ( iLoc > 0) /* ===== One isolated NonZero Value - no need to store a Window Len */ + if (iLoc > 0) /* ===== One isolated NonZero Value - no need to store a + Window Len */ { - vec[iLoc] = ((int *)raw)[0]; - raw += iShift; - } - else /* ===== New Window of size iLen */ + vec[iLoc] = ((int *)raw)[0]; + raw += iShift; + } else /* ===== New Window of size iLen */ { - int *adr = vec + (-iLoc); /* Start of the Windows in the output vector */ - iLen = *raw++; /* Length of the Window */ - - if (iLen > 0) /* ===== Non Constant Values */ - { - MEMCOPY( (int *)raw, adr, iLen, int); - - raw += iShift*iLen; - } - else /* ===== Constant Value : only one value is stored */ - { - iLen = - iLen; - int ValCst = *((int *)(raw)); raw += iShift; - do (*(adr++) = ValCst); while (--iLen > 0); - } + int *adr = + vec + (-iLoc); /* Start of the Windows in the output vector */ + iLen = *raw++; /* Length of the Window */ + + if (iLen > 0) /* ===== Non Constant Values */ + { + MEMCOPY((int *)raw, adr, iLen, int); + + raw += iShift * iLen; + } else /* ===== Constant Value : only one value is stored */ + { + iLen = -iLen; + int ValCst = *((int *)(raw)); + raw += iShift; + do + (*(adr++) = ValCst); + while (--iLen > 0); + } } - } while ( --NWin > 0); + } while (--NWin > 0); - return (char*)vec; + return (char *)vec; } +char *ReadWindowedSparseBufferShort(int *raw, int *size, short *vec) { -char* ReadWindowedSparseBufferShort(int *raw, int *size, short *vec){ - - int iShift = sizeof(short)/sizeof(int); + int iShift = sizeof(short) / sizeof(int); *size = *raw++; int NWin = *raw++; short *adr = vec; - MEM_ZERO( adr, *size*sizeof(short)); + MEM_ZERO(adr, *size * sizeof(short)); int iLoc, iLen; - if ( NWin > 0) - do { + if (NWin > 0) + do { - /* ===== We read the location of the new Window */ - iLoc = *raw++; /* ===== iLoc = Where start the next window */ + /* ===== We read the location of the new Window */ + iLoc = *raw++; /* ===== iLoc = Where start the next window */ - if ( iLoc > 0) /* ===== One isolated NonZero Value - no need to store a Window Len */ + if (iLoc > 0) /* ===== One isolated NonZero Value - no need to store a + Window Len */ { - vec[iLoc] = ((short *)raw)[0]; - raw += iShift; - } - else /* ===== New Window of size iLen */ + vec[iLoc] = ((short *)raw)[0]; + raw += iShift; + } else /* ===== New Window of size iLen */ { - short *adr = vec + (-iLoc); /* Start of the Windows in the output vector */ - iLen = *raw++; /* Length of the Window */ - - if (iLen > 0) /* ===== Non Constant Values */ - { - MEMCOPY( (short *)raw, adr, iLen, short); - - raw += iShift*iLen; - } - else /* ===== Constant Value : only one value is stored */ - { - iLen = - iLen; - short ValCst = *((short *)(raw)); raw += iShift; - do (*(adr++) = ValCst); while (--iLen > 0); - } + short *adr = + vec + (-iLoc); /* Start of the Windows in the output vector */ + iLen = *raw++; /* Length of the Window */ + + if (iLen > 0) /* ===== Non Constant Values */ + { + MEMCOPY((short *)raw, adr, iLen, short); + + raw += iShift * iLen; + } else /* ===== Constant Value : only one value is stored */ + { + iLen = -iLen; + short ValCst = *((short *)(raw)); + raw += iShift; + do + (*(adr++) = ValCst); + while (--iLen > 0); + } } - } while ( --NWin > 0); + } while (--NWin > 0); - return (char*)vec; + return (char *)vec; } - // read a record and return the pointer to the array -void* read_record(const char* filename, int64_t ptr, int* prec_flag, int* type_flag, - int* size, int* out_bufsize){ +void *read_record(const char *filename, int64_t ptr, int *prec_flag, + int *type_flag, int *size, int *out_bufsize) { int bsparse_flag, wsparse_flag, zlib_flag; - ifstream binFile (filename, ios::in | ios::binary); - binFile.seekg(ptr*4); - int bufsize = read_header(&binFile, &bsparse_flag, &wsparse_flag, - &zlib_flag, prec_flag, type_flag); + ifstream binFile(filename, ios::in | ios::binary); + binFile.seekg(ptr * 4); + int bufsize = read_header(&binFile, &bsparse_flag, &wsparse_flag, &zlib_flag, + prec_flag, type_flag); *size = bufsize; // always read record - char *raw = new char[4*bufsize]; - binFile.read(raw, 4*bufsize); - *out_bufsize = bufsize + 3; // include header and footer - - if (bsparse_flag){ - if (*type_flag){ - if (*prec_flag){ - raw = ReadShortBsparseRecord((int*)raw, size); - } else{ - raw = ReadBsparseRecord((int*)raw, size); + char *raw = new char[4 * bufsize]; + binFile.read(raw, 4 * bufsize); + *out_bufsize = bufsize + 3; // include header and footer + + if (bsparse_flag) { + if (*type_flag) { + if (*prec_flag) { + raw = ReadShortBsparseRecord((int *)raw, size); + } else { + raw = ReadBsparseRecord((int *)raw, size); } - } else{ // a float/double - if (*prec_flag){ - raw = ReadBsparseRecord((float*)raw, size); - } else{ - raw = ReadBsparseRecord((double*)raw, size); + } else { // a float/double + if (*prec_flag) { + raw = ReadBsparseRecord((float *)raw, size); + } else { + raw = ReadBsparseRecord((double *)raw, size); } } } else if (wsparse_flag) { - if (*type_flag){ - if (*prec_flag){ - raw = ReadWindowedSparseBuffer((short*)raw, size); - } else{ - raw = ReadWindowedSparseBuffer((int*)raw, size); + if (*type_flag) { + if (*prec_flag) { + raw = ReadWindowedSparseBuffer((short *)raw, size); + } else { + raw = ReadWindowedSparseBuffer((int *)raw, size); } - } else{ // a float/double - if (*prec_flag){ - raw = ReadWindowedSparseBuffer((float*)raw, size); - } else{ - raw = ReadWindowedSparseBuffer((double*)raw, size); + } else { // a float/double + if (*prec_flag) { + raw = ReadWindowedSparseBuffer((float *)raw, size); + } else { + raw = ReadWindowedSparseBuffer((double *)raw, size); } } - } return raw; } - - // populate arr with a record -// This function differs from read_record as it must be supplied with ``arr``, which must be sized properly to support the data coming from the file. -void read_record_stream(ifstream* file, int64_t loc, void* arr, int* prec_flag, - int* type_flag, int* size){ +// This function differs from read_record as it must be supplied with ``arr``, +// which must be sized properly to support the data coming from the file. +void read_record_stream(ifstream *file, int64_t loc, void *arr, int *prec_flag, + int *type_flag, int *size) { // seek to data location if supplied with a pointer - if (loc >= 0){ - file->seekg(loc*4); + if (loc >= 0) { + file->seekg(loc * 4); } int bsparse_flag, wsparse_flag, zlib_flag; - int bufsize = read_header(file, &bsparse_flag, &wsparse_flag, - &zlib_flag, prec_flag, type_flag); + int bufsize = read_header(file, &bsparse_flag, &wsparse_flag, &zlib_flag, + prec_flag, type_flag); *size = bufsize; // always read record - if (bufsize <= 0){ + if (bufsize <= 0) { return; } - char *raw = new char[4*bufsize]; + char *raw = new char[4 * bufsize]; - if (bsparse_flag){ + if (bsparse_flag) { // write to temporary record - file->read(raw, 4*bufsize); - - if (*type_flag){ - if (*prec_flag){ - ReadShortBsparseRecordToVec((int*)raw, size, (short*)arr); - } else{ - ReadBsparseRecordToVec((int*)raw, size, (int*)arr); + file->read(raw, 4 * bufsize); + + if (*type_flag) { + if (*prec_flag) { + ReadShortBsparseRecordToVec((int *)raw, size, (short *)arr); + } else { + ReadBsparseRecordToVec((int *)raw, size, (int *)arr); } - } else{ // a float or a double - if (*prec_flag){ - ReadBsparseRecordToVec((int*)raw, size, (float*)arr); - } else{ - ReadBsparseRecordToVec((int*)raw, size, (double*)arr); + } else { // a float or a double + if (*prec_flag) { + ReadBsparseRecordToVec((int *)raw, size, (float *)arr); + } else { + ReadBsparseRecordToVec((int *)raw, size, (double *)arr); } } } else if (wsparse_flag) { - file->read(raw, 4*bufsize); - if (*type_flag){ - if (*prec_flag){ - ReadWindowedSparseBufferShort((int*)raw, size, (short*)arr); - } else{ - ReadWindowedSparseBufferInt((int*)raw, size, (int*)arr); + file->read(raw, 4 * bufsize); + if (*type_flag) { + if (*prec_flag) { + ReadWindowedSparseBufferShort((int *)raw, size, (short *)arr); + } else { + ReadWindowedSparseBufferInt((int *)raw, size, (int *)arr); } - } else{ // a float/double - if (*prec_flag){ - ReadWindowedSparseBufferFloat((int*)raw, size, (float*)arr); - } else{ - ReadWindowedSparseBufferDouble((int*)raw, size, (double*)arr); + } else { // a float/double + if (*prec_flag) { + ReadWindowedSparseBufferFloat((int *)raw, size, (float *)arr); + } else { + ReadWindowedSparseBufferDouble((int *)raw, size, (double *)arr); } } - } else {// write directly to the array - file->read((char*)arr, 4*bufsize); - } + } else { // write directly to the array + file->read((char *)arr, 4 * bufsize); + } delete[] raw; - } - // read a record given a file stream -void* read_record_fid(fstream* file, int64_t loc, int* prec_flag, int* type_flag, - int* size, int* out_bufsize){ +void *read_record_fid(fstream *file, int64_t loc, int *prec_flag, + int *type_flag, int *size, int *out_bufsize) { int bsparse_flag, wsparse_flag, zlib_flag; // seek to data location if supplied with a pointer - file->seekg(loc*4); - int bufsize = read_header(file, &bsparse_flag, &wsparse_flag, - &zlib_flag, prec_flag, type_flag); + file->seekg(loc * 4); + int bufsize = read_header(file, &bsparse_flag, &wsparse_flag, &zlib_flag, + prec_flag, type_flag); *size = bufsize; // always read record - char *raw = new char[4*bufsize]; - file->read(raw, 4*bufsize); - *out_bufsize = bufsize + 3; // include header and footer - - if (bsparse_flag){ - if (*type_flag){ - if (*prec_flag){ - raw = ReadShortBsparseRecord((int*)raw, size); - } else{ - raw = ReadBsparseRecord((int*)raw, size); + char *raw = new char[4 * bufsize]; + file->read(raw, 4 * bufsize); + *out_bufsize = bufsize + 3; // include header and footer + + if (bsparse_flag) { + if (*type_flag) { + if (*prec_flag) { + raw = ReadShortBsparseRecord((int *)raw, size); + } else { + raw = ReadBsparseRecord((int *)raw, size); } - } else{ // a float/double - if (*prec_flag){ - raw = ReadBsparseRecord((float*)raw, size); - } else{ - raw = ReadBsparseRecord((double*)raw, size); + } else { // a float/double + if (*prec_flag) { + raw = ReadBsparseRecord((float *)raw, size); + } else { + raw = ReadBsparseRecord((double *)raw, size); } } } else if (wsparse_flag) { - if (*type_flag){ - if (*prec_flag){ - raw = ReadWindowedSparseBuffer((short*)raw, size); - } else{ - raw = ReadWindowedSparseBuffer((int*)raw, size); + if (*type_flag) { + if (*prec_flag) { + raw = ReadWindowedSparseBuffer((short *)raw, size); + } else { + raw = ReadWindowedSparseBuffer((int *)raw, size); } - } else{ // a float/double - if (*prec_flag){ - raw = ReadWindowedSparseBuffer((float*)raw, size); - } else{ - raw = ReadWindowedSparseBuffer((double*)raw, size); + } else { // a float/double + if (*prec_flag) { + raw = ReadWindowedSparseBuffer((float *)raw, size); + } else { + raw = ReadWindowedSparseBuffer((double *)raw, size); } } - } return raw; } - - -void read_nodes(const char* filename, int64_t ptrLOC, int nrec, int *nnum, - double *nodes){ +void read_nodes(const char *filename, int64_t ptrLOC, int nrec, int *nnum, + double *nodes) { // max buf size - char *raw = new char[68*4]; - ifstream binFile (filename, ios::in | ios::binary); - binFile.seekg(ptrLOC*4); + char *raw = new char[68 * 4]; + ifstream binFile(filename, ios::in | ios::binary); + binFile.seekg(ptrLOC * 4); // read remainder of buffer excluding initial bytes and last bytes int bufsize, n; int prec_flag, type_flag; - for (n=0; nseekg(ptr*4); - int bufsize = read_header(fs, &bsparse_flag, &wsparse_flag, - &zlib_flag, &prec_flag, &type_flag); + fs->seekg(ptr * 4); + int bufsize = read_header(fs, &bsparse_flag, &wsparse_flag, &zlib_flag, + &prec_flag, &type_flag); // need to perform size check // No compression allowed - if (bsparse_flag || wsparse_flag || zlib_flag){ + if (bsparse_flag || wsparse_flag || zlib_flag) { return 1; } @@ -651,63 +637,61 @@ int overwriteRecord(fstream* fs, int ptr, double* data){ // prec_flag 1 int32 or double // seek back to right after the header - fs->seekg(ptr*4 + 8); + fs->seekg(ptr * 4 + 8); - if (type_flag){ // int16 or int32 - return 1; // type not supported yet - } else { // either float or double - if (prec_flag){ // float + if (type_flag) { // int16 or int32 + return 1; // type not supported yet + } else { // either float or double + if (prec_flag) { // float // ideally we'd just copy this all at once // float* float_array[bufsize]; cout << "here" << endl; // std::copy(data, data + bufsize, *float_array); cout << "here" << endl; - // fs->write(reinterpret_cast(&float_array), bufsize*sizeof(float)); + // fs->write(reinterpret_cast(&float_array), + // bufsize*sizeof(float)); // however, this works float val; - for (int i=0; iwrite(reinterpret_cast(&val), sizeof(float)); } - } else { // must be double - fs->write(reinterpret_cast(&data), bufsize*sizeof(double)); + } else { // must be double + fs->write(reinterpret_cast(&data), + bufsize * sizeof(double)); } } return 0; - } - // overwrite an existing ansys record at location ptr -int overwriteRecordFloat(fstream* fs, int ptr, float* data){ +int overwriteRecordFloat(fstream *fs, int ptr, float *data) { int bsparse_flag, wsparse_flag, zlib_flag, prec_flag, type_flag; int size; // read the header // seek to data location if supplied with a pointer - fs->seekg(ptr*4); - int bufsize = read_header(fs, &bsparse_flag, &wsparse_flag, - &zlib_flag, &prec_flag, &type_flag); + fs->seekg(ptr * 4); + int bufsize = read_header(fs, &bsparse_flag, &wsparse_flag, &zlib_flag, + &prec_flag, &type_flag); // need to perform size check // No compression allowed - if (bsparse_flag || wsparse_flag || zlib_flag){ + if (bsparse_flag || wsparse_flag || zlib_flag) { return 1; } // seek back to right after the header - fs->seekg(ptr*4 + 8); - if (type_flag){ // either int16 or int32 - return 1; // type not supported yet - } else { // int16 or int32 - if (prec_flag){ // float - fs->write(reinterpret_cast(data), bufsize*sizeof(float)); - } else { // must be double - + fs->seekg(ptr * 4 + 8); + if (type_flag) { // either int16 or int32 + return 1; // type not supported yet + } else { // int16 or int32 + if (prec_flag) { // float + fs->write(reinterpret_cast(data), bufsize * sizeof(float)); + } else { // must be double } } return 0; - } diff --git a/ansys/mapdl/reader/cython/parsefull.c b/ansys/mapdl/reader/cython/parsefull.c index 78b5a667..ed4e60b7 100644 --- a/ansys/mapdl/reader/cython/parsefull.c +++ b/ansys/mapdl/reader/cython/parsefull.c @@ -8,316 +8,309 @@ // To add // - add big endianness compatibility +__inline int read_int(FILE *ptr) { + int n; + unsigned char c[4]; -__inline int read_int(FILE* ptr){ - int n; - unsigned char c[4]; + // Read and convert assuming small endian + fread(c, sizeof(int), 1, ptr); + n = (int)c[0] + ((int)c[1] << 8) + ((int)c[2] << 16) + ((int)c[3] << 24); - // Read and convert assuming small endian - fread(c, sizeof(int), 1, ptr); - n = (int)c[0] + ((int)c[1] << 8) + ((int)c[2]<<16) + ((int)c[3]<<24); - - return n; + return n; } -__inline double read_double(FILE* ptr){ - double val; - char str[8]; +__inline double read_double(FILE *ptr) { + double val; + char str[8]; - fread(str, sizeof(double), 1, ptr); - memcpy(&val, str, sizeof(double)); - return val; + fread(str, sizeof(double), 1, ptr); + memcpy(&val, str, sizeof(double)); + return val; } // returns an integer from a position along a char array -__inline int read_int_raw(char **c, int e){ - // Assuming small endian - int n; - memcpy(&n, (*c) + e, sizeof(int)); - return n; +__inline int read_int_raw(char **c, int e) { + // Assuming small endian + int n; + memcpy(&n, (*c) + e, sizeof(int)); + return n; } // returns a double from a position along a char array -__inline double read_double_raw(char **c, int e){ - double val; - memcpy(&val, (*c) + e, sizeof(double)); - return val; +__inline double read_double_raw(char **c, int e) { + double val; + memcpy(&val, (*c) + e, sizeof(double)); + return val; } - // index sorting int *array; -int cmp(const void *a, const void *b){ - int ia = *(int *)a; - int ib = *(int *)b; - return array[ia] < array[ib] ? -1 : array[ia] > array[ib]; +int cmp(const void *a, const void *b) { + int ia = *(int *)a; + int ib = *(int *)b; + return array[ia] < array[ib] ? -1 : array[ia] > array[ib]; } // Make index array -void make_index(int **idx, int size){ - int i; - *idx = (int*)malloc(size * sizeof(int)); +void make_index(int **idx, int size) { + int i; + *idx = (int *)malloc(size * sizeof(int)); - for(i=0; i 0){ - (*isfree)[i] = 1; - ++nfree; - } - else{ - (*isfree)[i] = 0; - } +int pop_isfree(FILE *ptr, int ptrDOF, int nNodes, int **isfree) { + int val, i; + int nfree = 0; + + for (i = 0; i < nNodes * 3; i += 1) { + val = read_int(ptr); + + if (val > 0) { + (*isfree)[i] = 1; + ++nfree; + } else { + (*isfree)[i] = 0; } + } - return nfree; + return nfree; } - //============================================================================= // Read an array from fortran file //============================================================================= int read_array(int **rows, int **cols, double **data, int *rref, int *cref, - int *isfree, int nterm, int neqn, FILE* ptr, int fileloc, - int *skipped){ - // Max array size - int i, j, val; - int e = 0; - int nitems, row, col; - double vald; - int c = 0; - int d = 0; - int inval = 0; - int val2, vald2; - - // Read entire matrix to memory - int nread = neqn*24 + nterm*12; - char *raw = (char*)malloc(nread*sizeof(char)); - fseek(ptr, (fileloc)*4, SEEK_SET); - fread(raw, sizeof(char), nread, ptr); - - for (i = 0; i //============================================================================= -// Fast string to integer convert to ANSYS formatted integers +// Fast string to integer convert to ANSYS formatted integers //============================================================================= -__inline int fast_atoi(char* raw, int intsz){ +__inline int fast_atoi(char *raw, int intsz) { int val; int c; val = 0; - for (c=0; c= '0' && *raw <= '9') { + } else if (*raw >= '0' && *raw <= '9') { val += (*raw++ - '0') * k; k *= 0.1; } @@ -95,81 +93,79 @@ __inline int ans_strtod(char *raw, int fltsz, double *arr){ // 1.0000000000000E-001 int evalue = 0; int esign = 1; - if (*raw == 'e' || *raw == 'E'){ + if (*raw == 'e' || *raw == 'E') { raw++; // skip "E" // always a sign of some sort - if (*raw == '-'){ + if (*raw == '-') { esign = -1; } - raw++; i++; i++; // skip E and sign + raw++; + i++; + i++; // skip E and sign /* printf(" %d<%d ", i, fltsz); */ - for (; i= '0' && *raw <= '9') { + } else if (*raw >= '0' && *raw <= '9') { val += (*raw++ - '0') * k; k *= 0.1; } @@ -179,33 +175,32 @@ static inline double ans_strtod2(char *raw, int fltsz){ // 1.0000000000000E-001 int evalue = 0; int esign = 1; - if (*raw == 'e' || *raw == 'E'){ + if (*raw == 'e' || *raw == 'E') { raw++; // skip "E" // always a sign of some sort - if (*raw == '-'){ + if (*raw == '-') { esign = -1; } - raw++; i++; i++; // skip E and sign - for (; i 20){ */ @@ -414,26 +412,29 @@ int read_eblock(char *raw, int *elem_off, int *elem, int nelem, int intsz, // Field 10: Not Used raw += intsz; - + // Field 11: Element number - elem[c++] = fast_atoi(raw, intsz); raw += intsz; + elem[c++] = fast_atoi(raw, intsz); + raw += intsz; /* printf("reading element %d\n", elem[c - 1]); */ // Need an additional value for consistency with other formats elem[c++] = 0; // Read nodes in element - for (j=0; j 10){ - for (j=nnode; j<20; j++){ + if (nnode < 20 && nnode > 10) { + for (j = nnode; j < 20; j++) { elem[c++] = 0; } } @@ -442,7 +443,6 @@ int read_eblock(char *raw, int *elem_off, int *elem, int nelem, int intsz, /* elem[c++] = 0; */ /* } */ /* } */ - } // update file position @@ -453,14 +453,13 @@ int read_eblock(char *raw, int *elem_off, int *elem, int nelem, int intsz, return c; } - // Simply write an array to disk as ASCII -int write_array_ascii(const char* filename, const double *arr, - const int nvalues){ - FILE * stream = fopen(filename, "w"); +int write_array_ascii(const char *filename, const double *arr, + const int nvalues) { + FILE *stream = fopen(filename, "w"); int i; - for (i=0; i -#include #include +#include #include +#include /* VTK numbering for vtk cells */ uint8_t VTK_EMPTY_CELL = 0; @@ -21,21 +21,19 @@ uint8_t VTK_QUADRATIC_PYRAMID = 27; uint8_t VTK_QUADRATIC_WEDGE = 26; uint8_t VTK_QUADRATIC_HEXAHEDRON = 25; - // Contains data for VTK UnstructuredGrid struct VtkData { int64_t *offset; int64_t *cells; uint8_t *celltypes; - int loc; // current position within cells - int *nref; // conversion between ansys and vtk numbering + int loc; // current position within cells + int *nref; // conversion between ansys and vtk numbering }; struct VtkData vtk_data; - // Populate offset, cell type, and prepare the cell array for a cell -static inline void add_cell(bool build_offset, int n_points, uint8_t celltype){ - if (build_offset){ +static inline void add_cell(bool build_offset, int n_points, uint8_t celltype) { + if (build_offset) { vtk_data.offset[0] = vtk_data.loc; vtk_data.offset++; } @@ -46,19 +44,17 @@ static inline void add_cell(bool build_offset, int n_points, uint8_t celltype){ return; } - - /* ============================================================================ * Store hexahedral element in vtk arrays. ANSYS elements are * ordered in the same manner as VTK. - * + * * VTK DOCUMENTATION * Linear Hexahedral * The hexahedron is defined by the eight points (0-7) where * (0,1,2,3) is the base of the hexahedron which, using the right * hand rule, forms a quadrilaterial whose normal points in the * direction of the opposite face (4,5,6,7). - * + * * Quadradic Hexahedral * The ordering of the twenty points defining the cell is point ids * (0-7, 8-19) where point ids 0-7 are the eight corner vertices of @@ -78,29 +74,30 @@ static inline void add_cell(bool build_offset, int n_points, uint8_t celltype){ * 18 (2, 6) * 19 (3, 7) */ -static inline void add_hex(bool build_offset, const int *elem, int nnode){ +static inline void add_hex(bool build_offset, const int *elem, int nnode) { int i; bool quad = nnode > 8; - if (quad){ + if (quad) { add_cell(build_offset, 20, VTK_QUADRATIC_HEXAHEDRON); } else { add_cell(build_offset, 8, VTK_HEXAHEDRON); } // Always add linear - for (i=0; i<8; i++){ + for (i = 0; i < 8; i++) { vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; } // translate connectivity - if (quad){ - for (i=8; i 8; - if (quad){ + if (quad) { add_cell(build_offset, 15, VTK_QUADRATIC_WEDGE); } else { add_cell(build_offset, 6, VTK_WEDGE); @@ -149,7 +144,7 @@ static inline void add_wedge(bool build_offset, const int *elem, int nnode){ vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[5]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; - if (quad){ // todo: check if index > nnode - 1 + if (quad) { // todo: check if index > nnode - 1 vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[9]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[8]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[11]]; @@ -164,56 +159,53 @@ static inline void add_wedge(bool build_offset, const int *elem, int nnode){ return; } - -static inline void add_pyr(bool build_offset, const int *elem, int nnode){ - int i; // counter +static inline void add_pyr(bool build_offset, const int *elem, int nnode) { + int i; // counter bool quad = nnode > 8; - if (quad){ + if (quad) { add_cell(build_offset, 13, VTK_QUADRATIC_PYRAMID); } else { add_cell(build_offset, 5, VTK_PYRAMID); } // [0, 1, 2, 3, 4, X, X, X] - for (i=0; i<5; i++){ - vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; + for (i = 0; i < 5; i++) { + vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; } - if (quad){ // todo: check if index > nnode - 1 - for (i=8; i<12; i++){ + if (quad) { // todo: check if index > nnode - 1 + for (i = 8; i < 12; i++) { vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; } - for (i=16; i<20; i++){ + for (i = 16; i < 20; i++) { vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; } } return; } - - /* ============================================================================ * Stores tetrahedral element in vtk arrays. ANSYS elements are * ordered the same as vtk elements. - * + * * VTK DOCUMENTATION: * Linear * The tetrahedron is defined by the four points (0-3); where (0,1,2) * is the base of the tetrahedron which, using the right hand rule, * forms a triangle whose normal points in the direction of the * fourth point. - * - * Quadradic + * + * Quadradic * The cell includes a mid-edge node on each of the size edges of the * tetrahedron. The ordering of the ten points defining the cell is * point ids (0-3,4-9) where ids 0-3 are the four tetra vertices; and * point ids 4-9 are the midedge nodes between: * (0,1), (1,2), (2,0), (0,3), (1,3), and (2,3) ============================================================================ */ -static inline void add_tet(bool build_offset, const int *elem, int nnode){ +static inline void add_tet(bool build_offset, const int *elem, int nnode) { bool quad = nnode > 8; - if (quad){ + if (quad) { add_cell(build_offset, 10, VTK_QUADRATIC_TETRA); } else { add_cell(build_offset, 4, VTK_TETRA); @@ -226,7 +218,7 @@ static inline void add_tet(bool build_offset, const int *elem, int nnode){ vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; - if (quad){ // todo: check if index > nnode - 1 + if (quad) { // todo: check if index > nnode - 1 vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[8]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[9]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[11]]; @@ -238,12 +230,11 @@ static inline void add_tet(bool build_offset, const int *elem, int nnode){ return; } - // ANSYS Tetrahedral style [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] -static inline void add_tet10(bool build_offset, const int *elem, int nnode){ - int i; // counter +static inline void add_tet10(bool build_offset, const int *elem, int nnode) { + int i; // counter bool quad = nnode > 4; - if (quad){ + if (quad) { add_cell(build_offset, 10, VTK_QUADRATIC_TETRA); } else { add_cell(build_offset, 4, VTK_TETRA); @@ -255,13 +246,13 @@ static inline void add_tet10(bool build_offset, const int *elem, int nnode){ vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[3]]; - if (quad){ - for (i=4; i ", i, elem[i]); */ /* printf("%i, \n", vtk_data.nref[elem[i]]); */ vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[i]]; @@ -300,8 +290,8 @@ static inline void add_quad(bool build_offset, const int *elem, bool is_quad){ return; } -void add_tri(bool build_offset, const int *elem, bool quad){ - if (quad){ +void add_tri(bool build_offset, const int *elem, bool quad) { + if (quad) { add_cell(build_offset, 6, VTK_QUADRATIC_TRIANGLE); } else { add_cell(build_offset, 3, VTK_TRIANGLE); @@ -312,7 +302,7 @@ void add_tri(bool build_offset, const int *elem, bool quad){ vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; - if (quad){ + if (quad) { vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[4]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[5]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[7]]; @@ -321,7 +311,8 @@ void add_tri(bool build_offset, const int *elem, bool quad){ return; } -void add_tri_missing_midside(bool build_offset, const int *elem, const int n_mid){ +void add_tri_missing_midside(bool build_offset, const int *elem, + const int n_mid) { add_cell(build_offset, 6, VTK_QUADRATIC_TRIANGLE); int i; @@ -336,23 +327,23 @@ void add_tri_missing_midside(bool build_offset, const int *elem, const int n_mid for (i = 0; i < n_mid && i < 3; i++) { vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[elem_indices[i]]]; } - for (i=3; i>n_mid; i--){ + for (i = 3; i > n_mid; i--) { vtk_data.cells[vtk_data.loc++] = -1; } return; } -static inline void add_line(bool build_offset, const int *elem, int nnode){ +static inline void add_line(bool build_offset, const int *elem, int nnode) { bool is_quad; - if (nnode > 2){ + if (nnode > 2) { is_quad = elem[2] > 0; } else { is_quad = false; } /* printf("is_quad, %i\n", is_quad); */ - if (is_quad){ + if (is_quad) { add_cell(build_offset, 3, VTK_QUADRATIC_EDGE); } else { add_cell(build_offset, 2, VTK_LINE); @@ -361,27 +352,24 @@ static inline void add_line(bool build_offset, const int *elem, int nnode){ // edge nodes vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[1]]; - if (is_quad){ + if (is_quad) { vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[2]]; } return; } - -static inline void add_point(bool build_offset, const int *elem){ +static inline void add_point(bool build_offset, const int *elem) { add_cell(build_offset, 1, VTK_VERTEX); vtk_data.cells[vtk_data.loc++] = vtk_data.nref[elem[0]]; return; } - - /* Stores wedge element in vtk arrays. ANSYS elements are ordered differently than vtk elements. ANSYS orders counter-clockwise and VTK orders clockwise - + VTK DOCUMENTATION Linear Wedge The wedge is defined by the six points (0-5) where (0,1,2) is the @@ -389,7 +377,7 @@ base of the wedge which, using the right hand rule, forms a triangle whose normal points outward (away from the triangular face (3,4,5)). -Quadradic Wedge +Quadradic Wedge The ordering of the fifteen points defining the cell is point ids (0-5,6-14) where point ids 0-5 are the six corner vertices of the wedge, defined analogously to the six @@ -397,20 +385,18 @@ points in vtkWedge (points (0,1,2) form the base of the wedge which, using the right hand rule, forms a triangle whose normal points away from the triangular face (3,4,5)); followed by nine midedge nodes (6-14). Note that these midedge nodes correspond lie -on the edges defined by : +on the edges defined by : (0,1), (1,2), (2,0), (3,4), (4,5), (5,3), (0,3), (1,4), (2,5) */ - - /* ============================================================================ * function: ans_to_vtk * Convert raw ANSYS elements to a VTK UnstructuredGrid format. - * + * * Parameters * ---------- * nelem : Number of elements. - * + * * elem: Array of elements * Each element contains 10 items plus the nodes belonging to the * element. The first 10 items are: @@ -440,35 +426,35 @@ on the edges defined by : * nnode : Number of nodes. * * nnum : ANSYS Node numbering - * + * * build_offset: Enable, disable populating offset array - * + * * Returns (Given as as parameters) * ------- * offset : VTK offset array to populate - * + * * cells : VTK cell connectivity - * + * * celltypes: VTK cell types - * + * * ==========================================================================*/ int ans_to_vtk(const int nelem, const int *elem, const int *elem_off, - const int *type_ref, const int nnode, const int* nnum, - int64_t *offset, int64_t *cells, uint8_t *celltypes, - const int build_offset){ + const int *type_ref, const int nnode, const int *nnum, + int64_t *offset, int64_t *cells, uint8_t *celltypes, + const int build_offset) { bool is_quad; - int i; // counter - int nnode_elem; // number of nodes belonging to the element - int off; // location of the element nodes - int etype; // ANSYS element type number + int i; // counter + int nnode_elem; // number of nodes belonging to the element + int off; // location of the element nodes + int etype; // ANSYS element type number // index ansys node number to VTK C based compatible indexing // max node number should be last node // Consider using a hash table here instead - int *nref = (int*) malloc((nnum[nnode - 1] + 1) * sizeof(int)); - nref[0] = -1; // for missing midside nodes ANSYS uses a node number of 0 - for (i=0; i Date: Tue, 6 Jun 2023 16:41:59 -0600 Subject: [PATCH 3/4] drop build for 3.7 --- .github/workflows/testing-and-deployment.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/testing-and-deployment.yml b/.github/workflows/testing-and-deployment.yml index a5796c99..ed76294c 100644 --- a/.github/workflows/testing-and-deployment.yml +++ b/.github/workflows/testing-and-deployment.yml @@ -79,7 +79,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - python-version: ['3.7', '3.8', '3.9', '3.10', '3.11'] + python-version: ['3.8', '3.9', '3.10', '3.11'] os: [ubuntu-latest, windows-latest] steps: @@ -156,7 +156,7 @@ jobs: - name: Build wheels uses: pypa/cibuildwheel@v2.13.0 with: - python-version: 3.7 - 3.11 + python-version: 3.8 - 3.11 - name: List generated wheels run: | From df8a4b84479af6f018798b2fb1b90f606b6c7baf Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Tue, 6 Jun 2023 16:44:02 -0600 Subject: [PATCH 4/4] generate release notes on release --- .github/workflows/testing-and-deployment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/testing-and-deployment.yml b/.github/workflows/testing-and-deployment.yml index ed76294c..0605c80e 100644 --- a/.github/workflows/testing-and-deployment.yml +++ b/.github/workflows/testing-and-deployment.yml @@ -274,6 +274,7 @@ jobs: - name: Release uses: softprops/action-gh-release@v1 with: + generate_release_notes: true files: | ./**/*.whl ./**/*.zip