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