1 /*
2 * $Id: yhex.c,v 1.1 2005-09-18 22:05:51 dhmunro Exp $
3 * Yorick-specific interface routines, see hex.i
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 #include "hex.h"
12 #include "regul.h"
13 #include "ydata.h"
14 #include "yio.h"
15 #include "pstdlib.h"
16
17 extern BuiltIn Y_hex_mesh, Y_hex5_track, Y_hex24f_track, Y_hex24b_track;
18 extern BuiltIn Y_reg_track, Y_hex_query;
19
20 /* Implement HX_mesh as a foreign Yorick data type. */
21 typedef struct YHX_mesh YHX_mesh;
22 struct YHX_mesh {
23 int references; /* reference counter */
24 Operations *ops; /* virtual function table */
25 HX_mesh mesh;
26 TK_result *result;
27 };
28
29 extern YHX_mesh *new_YHX(double *xyz, long *bound, long nbnds, HX_blkbnd *bnds,
30 long nblks, HX_block *blks, long start);
31 extern void free_YHX(void *yhx); /* ******* Use Unref(dm) ******* */
32 extern YHX_mesh *YGet_YHX_mesh(Symbol *s);
33 extern Operations yhx_mesh_ops;
34
Y_hex_mesh(int nArgs)35 void Y_hex_mesh(int nArgs)
36 {
37 Symbol *s= sp-nArgs+1;
38 if (nArgs==7) {
39 double *xyz= YGet_D(s, 0, (Dimension **)0);
40 long *bound= YGet_L(s+1, 0, (Dimension **)0);
41 long nbnds= YGetInteger(s+2);
42 void **pbnds= YGet_P(s+3, 1, (Dimension **)0);
43 long nblks= YGetInteger(s+4);
44 void **pblks= YGet_P(s+5, 0, (Dimension **)0);
45 long start= YGetInteger(s+6);
46 if (!pblks) YError("hex_mesh blks parameter bad");
47 PushDataBlock(new_YHX(xyz, bound, pbnds? nbnds : 0,
48 pbnds? *pbnds : 0, nblks, *pblks, start));
49 } else {
50 YError("hex_mesh needs exactly seven arguments");
51 }
52 }
53
Y_hex_query(int nArgs)54 void Y_hex_query(int nArgs)
55 {
56 Symbol *s= sp-nArgs+1;
57 YHX_mesh *mesh;
58 if (nArgs<1 || nArgs>5) YError("hex_query needs 1-5 arguments");
59 if (s->ops==&referenceSym) ReplaceRef(s);
60 if (s->ops!=&dataBlockSym || s->value.db->ops!=&yhx_mesh_ops)
61 YError("hex_query 1st argument must be a hex mesh");
62 mesh= (YHX_mesh *)s->value.db;
63 if (++s<=sp) {
64 Symbol sy;
65 long index= YGet_Ref(s);
66 sy.ops= &dataBlockSym;
67 sy.value.db= Pointee(mesh->mesh.xyz);
68 YPut_Result(&sy, index);
69 if (++s<=sp) {
70 index= YGet_Ref(s);
71 sy.value.db= Pointee(mesh->mesh.bound);
72 YPut_Result(&sy, index);
73 if (++s<=sp) {
74 index= YGet_Ref(s);
75 sy.value.db= Pointee(mesh->mesh.bnds);
76 YPut_Result(&sy, index);
77 if (++s<=sp) {
78 index= YGet_Ref(s);
79 sy.value.db= Pointee(mesh->mesh.blks);
80 YPut_Result(&sy, index);
81 }
82 }
83 }
84 }
85 PushLongValue(mesh->mesh.start);
86 }
87
new_YHX(double * xyz,long * bound,long nbnds,HX_blkbnd * bnds,long nblks,HX_block * blks,long start)88 YHX_mesh *new_YHX(double *xyz, long *bound, long nbnds, HX_blkbnd *bnds,
89 long nblks, HX_block *blks, long start)
90 {
91 YHX_mesh *mesh= p_malloc(sizeof(YHX_mesh));
92 Array *ary;
93 mesh->references= 0;
94 mesh->ops= &yhx_mesh_ops;
95 mesh->result= 0;
96 mesh->mesh.xyz= (void *)xyz;
97 mesh->mesh.orient= 0;
98 mesh->mesh.stride= 0;
99 mesh->mesh.bound= (void *)bound;
100 mesh->mesh.nbnds= nbnds;
101 mesh->mesh.bnds= bnds;
102 mesh->mesh.nblks= nblks;
103 mesh->mesh.blks= blks;
104 mesh->mesh.block= 0;
105 mesh->mesh.start= start;
106 if (xyz) {
107 ary= Pointee(xyz);
108 Ref(ary);
109 }
110 if (bound) {
111 ary= Pointee(bound);
112 Ref(ary);
113 }
114 if (bnds) {
115 ary= Pointee(bnds);
116 Ref(ary);
117 }
118 if (blks) {
119 ary= Pointee(blks);
120 Ref(ary);
121 }
122 return mesh;
123 }
124
free_YHX(void * yhx)125 void free_YHX(void *yhx)
126 {
127 YHX_mesh *mesh= yhx;
128 Array *ary;
129 TK_result *result= mesh->result;
130 mesh->result= 0;
131 if (result) ray_free(result);
132 ary= mesh->mesh.xyz? Pointee(mesh->mesh.xyz) : 0;
133 mesh->mesh.xyz= 0;
134 if (ary) Unref(ary);
135 ary= mesh->mesh.bound? Pointee(mesh->mesh.bound) : 0;
136 mesh->mesh.bound= 0;
137 if (ary) Unref(ary);
138 ary= mesh->mesh.bnds? Pointee(mesh->mesh.bnds) : 0;
139 mesh->mesh.bnds= 0;
140 if (ary) Unref(ary);
141 ary= mesh->mesh.blks? Pointee(mesh->mesh.blks) : 0;
142 mesh->mesh.blks= 0;
143 if (ary) Unref(ary);
144 p_free(mesh);
145 }
146
147 static void hex_tracker(int nArgs, int which);
148
Y_hex5_track(int nArgs)149 void Y_hex5_track(int nArgs)
150 {
151 hex_tracker(nArgs, 0);
152 }
153
Y_hex24f_track(int nArgs)154 void Y_hex24f_track(int nArgs)
155 {
156 hex_tracker(nArgs, 1);
157 }
158
Y_hex24b_track(int nArgs)159 void Y_hex24b_track(int nArgs)
160 {
161 hex_tracker(nArgs, 2);
162 }
163
164 static double *normalize_rays(double **p, long n);
165
hex_tracker(int nArgs,int which)166 static void hex_tracker(int nArgs, int which)
167 {
168 YHX_mesh *mesh;
169 double *p, *q;
170 Dimension *dims;
171 long index, dlist[10], n;
172 int ndims, i;
173 TK_result *result;
174 Array *c, *s;
175 if (nArgs!=3) YError("hexN_track takes exactly 3 arguments");
176 mesh= YGet_YHX_mesh(sp-2);
177 p= YGet_D(sp-1, 0, &dims);
178 index= YGet_Ref(sp);
179 Drop(1); /* s argument is just an output */
180 ndims= YGet_dims(dims, dlist, 10);
181 if (ndims>10 || ndims<2 || dlist[0]!=3 || dlist[ndims-1]!=2)
182 YError("hexN_track rays must be 3 x ray_dims x 2 array of [p,q]");
183 for (n=1,i=1 ; i<ndims-1 ; i++) n*= dlist[i];
184 q= normalize_rays(&p, n);
185 result= mesh->result;
186 if (result) ray_reset(result);
187 else result= mesh->result= ray_result();
188 if (!which)
189 hex5_rays(&mesh->mesh, n, (void *)p, (void *)q, result);
190 else if (which==1)
191 hex24_rays(&mesh->mesh, n, (void *)p, (void *)q, 0, result);
192 else
193 hex24_rays(&mesh->mesh, n, (void *)p, (void *)q, 1, result);
194 n= ray_collect(result, (long *)0, (double *)0, 1L);
195 dims= tmpDims;
196 tmpDims= 0;
197 FreeDimension(dims);
198 tmpDims= NewDimension(n, 1L, (Dimension *)0);
199 s= PushDataBlock(NewArray(&doubleStruct, tmpDims));
200 YPut_Result(sp, index);
201 c= PushDataBlock(NewArray(&longStruct, tmpDims));
202 ray_collect(result, c->value.l, s->value.d, 1L);
203 mesh->result= 0;
204 ray_free(result);
205 }
206
Y_reg_track(int nArgs)207 void Y_reg_track(int nArgs)
208 {
209 YHX_mesh *mesh;
210 double *p, *q, *xyz[3];
211 Dimension *dims;
212 long index, dlist[10], n, nxyz[3];
213 int ndims, i;
214 TK_result *result;
215 Array *c, *s;
216 if (nArgs!=5) YError("reg_track takes exactly 5 arguments");
217 for (i=0 ; i<3 ; i++) {
218 xyz[i]= YGet_D(sp-4+i, 0, &dims);
219 if (YGet_dims(dims, dlist, 2)!=1 || dlist[0]<2)
220 YError("reg_track x,y,z arguments must be 1D with >=2 elements");
221 nxyz[i]= dlist[0];
222 }
223 p= YGet_D(sp-1, 0, &dims);
224 index= YGet_Ref(sp);
225 Drop(1); /* s argument is just an output */
226 ndims= YGet_dims(dims, dlist, 10);
227 if (ndims>10 || ndims<2 || dlist[0]!=3 || dlist[ndims-1]!=2)
228 YError("reg_track rays must be 3 x ray_dims x 2 array of [p,q]");
229 for (n=1,i=1 ; i<ndims-1 ; i++) n*= dlist[i];
230 q= normalize_rays(&p, n);
231 /* push dummy mesh onto stack to hold result, allowing proper
232 * deallocation in case of an interrupt */
233 mesh= PushDataBlock(new_YHX((double *)0, (long *)0, 0, (HX_blkbnd *)0,
234 0, (HX_block *)0, 0));
235 result= mesh->result= ray_result();
236 reg_rays(nxyz, xyz, n, (void *)p, (void *)q, result);
237 n= ray_collect(result, (long *)0, (double *)0, 1L);
238 dims= tmpDims;
239 tmpDims= 0;
240 FreeDimension(dims);
241 tmpDims= NewDimension(n, 1L, (Dimension *)0);
242 s= PushDataBlock(NewArray(&doubleStruct, tmpDims));
243 YPut_Result(sp, index);
244 Drop(1);
245 c= PushDataBlock(NewArray(&longStruct, tmpDims));
246 ray_collect(result, c->value.l, s->value.d, 1L);
247 }
248
normalize_rays(double ** p,long n)249 static double *normalize_rays(double **p, long n)
250 {
251 extern double sqrt(double);
252 double *q, qnrm, qtst;
253 long i;
254 Array *array= (Array *)sp->value.db;
255 if (sp->ops!=&dataBlockSym || !array->ops->isArray)
256 YError("(BUG) normalize_rays failed");
257 if (array->references) {
258 /* copy non-temporary arrays to avoid unexpected aliasing */
259 Array *result= PushDataBlock(NewArray(array->type.base,
260 array->type.dims));
261 array->type.base->Copy(array->type.base, result->value.c,
262 array->value.c, array->type.number);
263 PopTo(sp-1);
264 *p= result->value.d;
265 }
266 q= *p + 3*n;
267 for (i=0 ; i<3*n ; i+=3) {
268 qnrm= q[i]<0.? -q[i] : q[i];
269 qtst= q[i+1]<0.? -q[i+1] : q[i+1];
270 if (qtst > qnrm) qnrm= qtst;
271 qtst= q[i+2]<0.? -q[i+2] : q[i+2];
272 if (qtst > qnrm) qnrm= qtst;
273 if (qnrm) {
274 qnrm= 1.0/qnrm;
275 q[i]*= qnrm;
276 q[i+1]*= qnrm;
277 q[i+2]*= qnrm;
278 qnrm= 1.0/sqrt(q[i]*q[i]+q[i+1]*q[i+1]+q[i+2]*q[i+2]);
279 q[i]*= qnrm;
280 q[i+1]*= qnrm;
281 q[i+2]*= qnrm;
282 } else {
283 q[i]= q[i+1]= 0.0;
284 q[i+2]= 1.0;
285 }
286 }
287 return q;
288 }
289
290 /*--------------------------------------------------------------------------*/
291 /* Define opaque mesh object YHX_mesh as a foreign Yorick data type. */
292
293 extern PromoteOp PromXX;
294 extern UnaryOp ToAnyX, NegateX, ComplementX, NotX, TrueX;
295 extern BinaryOp AddX, SubtractX, MultiplyX, DivideX, ModuloX, PowerX;
296 extern BinaryOp EqualX, NotEqualX, GreaterX, GreaterEQX;
297 extern BinaryOp ShiftLX, ShiftRX, OrX, AndX, XorX;
298 extern BinaryOp AssignX, MatMultX;
299 extern UnaryOp EvalX, SetupX, PrintX;
300 extern MemberOp GetMemberX;
301
302 static UnaryOp PrintYHX;
303
304 Operations yhx_mesh_ops = {
305 &free_YHX, T_OPAQUE, 0, T_STRING, "Hex-Mesh",
306 {&PromXX, &PromXX, &PromXX, &PromXX, &PromXX, &PromXX, &PromXX, &PromXX},
307 &ToAnyX, &ToAnyX, &ToAnyX, &ToAnyX, &ToAnyX, &ToAnyX, &ToAnyX,
308 &NegateX, &ComplementX, &NotX, &TrueX,
309 &AddX, &SubtractX, &MultiplyX, &DivideX, &ModuloX, &PowerX,
310 &EqualX, &NotEqualX, &GreaterX, &GreaterEQX,
311 &ShiftLX, &ShiftRX, &OrX, &AndX, &XorX,
312 &AssignX, &EvalX, &SetupX, &GetMemberX, &MatMultX, &PrintYHX
313 };
314
PrintYHX(Operand * op)315 static void PrintYHX(Operand *op)
316 {
317 YHX_mesh *yhx= op->value;
318 char line[96];
319 ForceNewline();
320 sprintf(line, "hex mesh: %ld blocks, %ld nodes", yhx->mesh.nblks,
321 yhx->mesh.blks[yhx->mesh.nblks-1].final);
322 PrintFunc(line);
323 ForceNewline();
324 }
325
YGet_YHX_mesh(Symbol * s)326 YHX_mesh *YGet_YHX_mesh(Symbol *s)
327 {
328 if (s->ops==&referenceSym) ReplaceRef(s);
329 if (s->ops!=&dataBlockSym || s->value.db->ops!=&yhx_mesh_ops)
330 YError("expecting Hex-Mesh argument");
331 return (YHX_mesh *)s->value.db;
332 }
333
334 /*--------------------------------------------------------------------------*/
335