source: issm/trunk/src/m/extrusion/project3d.py@ 26744

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

merged trunk-jpl and trunk for revision 26742
File size: 8.3 KB
Line 
1import numpy as np
2from pairoptions import pairoptions
3
4
5def project3d(md, *args):
6 '''
7 PROJECT3D - vertically project a vector from 2d mesh
8
9 vertically project a vector from 2d mesh (split in noncoll and coll
10 areas) into a 3d mesh.
11 This vector can be a node vector of size (md.mesh.numberofvertices2d,
12 N/A) or an element vector of size (md.mesh.numberofelements2d, N/A).
13
14 arguments:
15 'vector': 2d vector
16 'type': 'element' or 'node' or 'poly'
17
18 options:
19 'layer' a layer number where vector should keep its values. If
20 not specified, all layers adopt the value of the 2d
21 vector.
22 'padding': default to 0 (value adopted by other 3d layers not
23 being projected.
24 'degree': degree of polynomials when extrude from bottom to the top
25
26 Examples:
27 extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN)
28 extruded_vector = project3d(md, 'vector', vector2d, 'type', 'element', 'padding', 0)
29 extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node')
30 '''
31
32 #some regular checks
33 if not md:
34 raise TypeError("bad usage")
35 if md.mesh.domaintype().lower() != '3d':
36 raise TypeError("input model is not 3d")
37
38 #retrieve parameters from options.
39 options = pairoptions(*args)
40 vector2d = options.getfieldvalue('vector') #mandatory
41 vectype = options.getfieldvalue('type') #mandatory
42 layer = options.getfieldvalue('layer', 0) #optional (do all layers otherwise)
43 paddingvalue = options.getfieldvalue('padding', 0) #0 by default
44 polyexponent = options.getfieldvalue('degree', 0) #0 by default, 0-degree polynomial
45
46 #Handle special case where vector2d is single element (differs from representation in MATLAB)
47 if isinstance(vector2d, (bool, int, float)):
48 projected_vector = vector2d
49
50 if np.size(vector2d) == 1:
51 projected_vector = vector2d
52
53 elif vectype.lower() == 'node':
54 #Initialize 3d vector
55 if np.ndim(vector2d) == 1:
56 if vector2d.shape[0] == md.mesh.numberofvertices2d:
57 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices))).astype(vector2d.dtype)
58 elif vector2d.shape[0] == md.mesh.numberofvertices2d + 1:
59 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices + 1))).astype(vector2d.dtype)
60 projected_vector[-1] = vector2d[-1]
61 vector2d = vector2d[:-1]
62 else:
63 raise TypeError("vector length not supported")
64 #Fill in
65 if layer == 0:
66 for i in range(md.mesh.numberoflayers):
67 projected_vector[(i * md.mesh.numberofvertices2d):((i + 1) * md.mesh.numberofvertices2d)] = vector2d
68 else:
69 projected_vector[((layer - 1) * md.mesh.numberofvertices2d):(layer * md.mesh.numberofvertices2d)] = vector2d
70 else:
71 if vector2d.shape[0] == md.mesh.numberofvertices2d:
72 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
73 elif vector2d.shape[0] == md.mesh.numberofvertices2d + 1:
74 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices + 1, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
75 projected_vector[-1, :] = vector2d[-1, :]
76 vector2d = vector2d[:-1, :]
77 else:
78 raise TypeError("vector length not supported")
79 #Fill in
80 if layer == 0:
81 for i in range(md.mesh.numberoflayers):
82 projected_vector[(i * md.mesh.numberofvertices2d):((i + 1) * md.mesh.numberofvertices2d), :] = vector2d
83 else:
84 projected_vector[((layer - 1) * md.mesh.numberofvertices2d):(layer * md.mesh.numberofvertices2d), :] = vector2d
85
86 elif vectype.lower() == 'element':
87 #Initialize 3d vector
88 if np.ndim(vector2d) == 1:
89 if vector2d.shape[0] == md.mesh.numberofelements2d:
90 projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements))).astype(vector2d.dtype)
91 elif vector2d.shape[0] == md.mesh.numberofelements2d + 1:
92 projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1))).astype(vector2d.dtype)
93 projected_vector[-1] = vector2d[-1]
94 vector2d = vector2d[:-1]
95 else:
96 raise TypeError("vector length not supported")
97 #Fill in
98 if layer == 0:
99 for i in range(md.mesh.numberoflayers - 1):
100 projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d)] = vector2d
101 else:
102 projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d)] = vector2d
103 else:
104 if vector2d.shape[0] == md.mesh.numberofelements2d:
105 projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
106 elif vector2d.shape[0] == md.mesh.numberofelements2d + 1:
107 projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
108 projected_vector[-1, :] = vector2d[-1, :]
109 vector2d = vector2d[:-1, :]
110 else:
111 raise TypeError("vector length not supported")
112 #Fill in
113 if layer == 0:
114 for i in range(md.mesh.numberoflayers - 1):
115 projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d), :] = vector2d
116 else:
117 projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d), :] = vector2d
118 elif vectype.lower() == 'poly':
119 #Initialize 3d vector
120 if np.ndim(vector2d) == 1:
121 if vector2d.shape[0] == md.mesh.numberofvertices2d:
122 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices))).astype(vector2d.dtype)
123 elif vector2d.shape[0] == md.mesh.numberofvertices2d + 1:
124 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices + 1))).astype(vector2d.dtype)
125 projected_vector[-1] = vector2d[-1]
126 vector2d = vector2d[:-1]
127 else:
128 raise TypeError("vector length not supported")
129 #Fill in
130 if layer == 0:
131 for i in range(md.mesh.numberoflayers - 1):
132 projected_vector[(i * md.mesh.numberofvertices2d):((i + 1) * md.mesh.numberofvertices2d)] = vector2d * (1.0 - (1.0 - i / (md.mesh.numberoflayers - 1.0))**polyexponent)
133 else:
134 projected_vector[((layer - 1) * md.mesh.numberofvertices2d):(layer * md.mesh.numberofvertices2d)] = vector2d * (1.0 - (1.0 - layer / (md.mesh.numberoflayers - 1.0))**polyexponent)
135 else:
136 if vector2d.shape[0] == md.mesh.numberofvertices2d:
137 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
138 elif vector2d.shape[0] == md.mesh.numberofvertices2d + 1:
139 projected_vector = (paddingvalue * np.ones((md.mesh.numberofvertices + 1, np.size(vector2d, axis=1)))).astype(vector2d.dtype)
140 projected_vector[-1, :] = vector2d[-1, :]
141 vector2d = vector2d[:-1, :]
142 else:
143 raise TypeError("vector length not supported")
144 #Fill in
145 if layer == 0:
146 for i in range(md.mesh.numberoflayers - 1):
147 projected_vector[(i * md.mesh.numberofvertices2d):((i + 1) * md.mesh.numberofvertices2d), :] = vector2d * (1.0 - (1.0 - i / (md.mesh.numberoflayers - 1.0))**polyexponent)
148 else:
149 projected_vector[((layer - 1) * md.mesh.numberofvertices2d):(layer * md.mesh.numberofvertices2d), :] = vector2d * (1.0 - (1.0 - layer / (md.mesh.numberoflayers - 1.0))**polyexponent)
150 else:
151 raise TypeError("project3d error message: unknown projection type")
152
153 return projected_vector
Note: See TracBrowser for help on using the repository browser.