-
Notifications
You must be signed in to change notification settings - Fork 89
/
mesh_factory_ext.pyx
89 lines (70 loc) · 2.29 KB
/
mesh_factory_ext.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#cython: wraparound=False, boundscheck=False, cdivision=True, profile=False, nonecheck=False, overflowcheck=False, cdivision_warnings=False, unraisable_tracebacks=False
import cython
# import both numpy and the Cython declarations for numpy
import numpy as np
cimport numpy as np
def rectangular_cross_construct(np.ndarray[double, ndim=1, mode="c"] params not None,\
np.ndarray[double, ndim=1, mode="c"] origin not None,\
np.ndarray[double, ndim=2, mode="c"] points not None,\
np.ndarray[long, ndim=2, mode="c"] elements not None):
cdef int m, n, i, j, v1, v2 ,v3 ,v4, v5
cdef int numPoints, numElements
cdef double len1, len2, delta1, delta2, x, y
m = int(params[0])
n = int(params[1])
len1 = params[2]
len2 = params[3]
cdef np.ndarray[long, ndim=2, mode="c"] vertices = np.ascontiguousarray(np.zeros((m+1,n+1),dtype=np.int))
delta1 = len1/m
delta2 = len2/n
numPoints = 0
for i in xrange(m+1):
for j in xrange(n+1):
vertices[i,j] = numPoints
points[numPoints,0] = i*delta1 + origin[0]
points[numPoints,1] = j*delta2 + origin[1]
numPoints += 1
boundary = {}
numElements = 0
for i in xrange(m):
for j in xrange(n):
v1 = vertices[i,j+1]
v2 = vertices[i,j]
v3 = vertices[i+1,j+1]
v4 = vertices[i+1,j]
x = (points[v1,0] + points[v2,0] + points[v3,0] + points[v4,0])*0.25
y = (points[v1,1] + points[v2,1] + points[v3,1] + points[v4,1])*0.25
# Create centre point
v5 = numPoints
points[numPoints,0] = x
points[numPoints,1] = y
numPoints += 1
# Create left triangle
if i == 0:
boundary[(numElements,1)] = "left"
elements[numElements,0] = v2
elements[numElements,1] = v5
elements[numElements,2] = v1
numElements += 1
# Create bottom triangle
if j == 0:
boundary[(numElements,1)] = "bottom"
elements[numElements,0] = v4
elements[numElements,1] = v5
elements[numElements,2] = v2
numElements += 1
# Create right triangle
if i == m-1:
boundary[(numElements,1)] = "right"
elements[numElements,0] = v3
elements[numElements,1] = v5
elements[numElements,2] = v4
numElements += 1
# Create top triangle
if j == n-1:
boundary[(numElements,1)] = "top"
elements[numElements,0] = v1
elements[numElements,1] = v5
elements[numElements,2] = v3
numElements += 1
return boundary