Revision Log

contmex5.c

Revision 495 (checked in by cbc, 12 years ago)	Initial import of Stark code.

Line
1	#include <math.h>
2	#include <stdio.h>
3	#include "mex.h"
4	#include "opnml_mex5_allocs.c"
5
6	/* PROTOTYPES */
7	void isopts(int,
8	double *,
9	double *,
10	int *,
11	double *,
12	double,
13	double **,
14	int *);
15
16	/************************************************************
17
18	#### ## ##### ###### # # ## # #
19	# # # # # # # # # # # #
20	# # # # ##### # # # # #
21	# ### ###### # # # ## # ###### #
22	# # # # # # ## ## # # #
23	#### # # # ###### # # # # #
24
25	************************************************************/
26
27	void mexFunction(int nlhs,
28	mxArray *plhs[],
29	int nrhs,
30	const mxArray *prhs[])
31	{
32
33	/* ---- contourmex will be called as :
34	cmat=contourmex(x,y,ele,q,cval); ---------------------------- */
35
36	int cnt,*ele,i,j,nn,ne;
37	double x, y, *q;
38	double cval,dele;
39	double *cmat,newcmat;
40	int nrl,nrh,ncl,nch ;
41
42	/* ---- check I/O arguments ----------------------------------------- */
43	if (nrhs != 5)
44	mexErrMsgTxt("contmex5 requires 5 input arguments.");
45	else if (nlhs != 1)
46	mexErrMsgTxt("contmex5 requires 1 output arguments.");
47
48	/* ---- dereference input arrays ------------------------------------ */
49	x=mxGetPr(prhs[0]);
50	y=mxGetPr(prhs[1]);
51	dele=mxGetPr(prhs[2]);
52	q=mxGetPr(prhs[3]);
53	cval=mxGetPr(prhs[4]);
54	nn=mxGetM(prhs[0]);
55	ne=mxGetM(prhs[2]);
56
57	/* ---- allocate space for int representation of dele &
58	convert double element representation to int &
59	shift node numbers toward 0 by 1 for proper indexing -------- */
60	ele=(int )mxIvector(0,3ne);
61	for (i=0;i<3*ne;i++)
62	ele[i]=((int)dele[i])-1;
63
64	/* ---- allocate space for contour list cmat following NRC
65	allocation style
66	double **mxDmatrix(int nrl,int nrh,int ncl,int nch)
67	cmat= mxDmatrix( 0, ne-1, 0, 3) ------------- */
68	nrl=0;
69	nrh=6*ne;
70	ncl=0;
71	nch=3;
72	cmat=(double **) mxDmatrix(nrl,nrh,ncl,nch);
73
74	isopts(ne,x,y,ele,q,cval[0],cmat,&cnt);
75
76	if(cnt!=0){
77	newcmat=(double ) mxDvector(0,4cnt-1);
78	for (j=0;j<4;j++)
79	for (i=0;i<cnt;i++)
80	newcmat[cnt*j+i]=cmat[i][j];
81	/* ---- Set elements of return matrix, pointed to by plhs[0] -------- */
82	plhs[0]=mxCreateDoubleMatrix(cnt,4,mxREAL);
83	mxSetPr(plhs[0],newcmat);
84	}
85	else {
86	plhs[0]=mxCreateDoubleMatrix(1,1,mxREAL);
87	mxSetPr(plhs[0],NULL);
88	}
89
90	/* ---- No need to free memory allocated with "mxCalloc"; MATLAB
91	does this automatically. The CMEX allocation functions in
92	"opnml_allocs.c" use mxCalloc. ----------------------------------- */
93	return;
94	}
95
96	/*----------------------------------------------------------------------
97
98	# #### #### ##### ##### ####
99	# # # # # # # #
100	# #### # # # # # ####
101	# # # # ##### # #
102	# # # # # # # # #
103	# #### #### # # ####
104
105	----------------------------------------------------------------------*/
106
107	#ifdef __STDC__
108	void isopts(int ne,
109	double x, double y,
110	int ele,double q,
111	double cval,
112	double *cmat,int cnt)
113	#else
114	void isopts(ne,x,y,ele,q,cval,cmat,cnt)
115	int ne;
116	double x, y, q,cval,*cmat;
117	int ele,cnt;
118	#endif
119	#define ELE(i,j,m) ele[i+m*j]
120	#define TOL 1.e-10
121	{
122
123	/* ---- ELE is defined to perform the following array
124	element extraction ELE(i,j,m) ele[i+mj] ---------------- /
125	int count=0,k;
126	int n0,n1,n2;
127	double s0,s1,s2;
128	double xa,xb,ya,yb;
129	double fac;
130	double x0,x1,x2;
131	double y0,y1,y2;
132
133	/* ---- Element loop to determine contours -------------------------- */
134	for(k=0;k<ne;k++){
135
136	/* ---- First arrange element k nodes such that n0,n1,n2 => s0<s1<s2 -*/
137	n0=ELE(k,0,ne);
138	s0=q[n0];
139	n1=ELE(k,1,ne);
140	s1=q[n1];
141	if (s1<s0){
142	n0=n1;
143	s0=s1;
144	n1=ELE(k,0,ne);
145	s1=q[n1];
146	}
147	n2=ELE(k,2,ne);
148	s2=q[n2];
149	if (s2<s1){
150	n2=n1;
151	s2=s1;
152	n1=ELE(k,2,ne);
153	s1=q[n1];
154	if(s1<s0){
155	n1=n0;
156	s1=s0;
157	n0=ELE(k,2,ne);
158	s0=q[n0];
159	}
160	}
161	x0=x[n0];
162	y0=y[n0];
163	x1=x[n1];
164	y1=y[n1];
165	x2=x[n2];
166	y2=y[n2];
167
168	if(cval<s0)goto L10; /* cval < element min ------------------ */
169	if(cval>s2)goto L10; /* cval > element max ------------------ */
170
171	if (fabs(s0-s1)<TOL&& /* Contour on side n0 -> n1 --- */
172	fabs(s0-cval)<TOL){
173	cmat[count][0]=x0;
174	cmat[count][1]=y0;
175	cmat[count][2]=x1;
176	cmat[count][3]=y1;
177	count++;
178	}
179	else if (fabs(s1-s2)<TOL&& /* Contour on side n1 -> n2 --- */
180	fabs(s1-cval)<TOL){
181	cmat[count][0]=x1;
182	cmat[count][1]=y1;
183	cmat[count][2]=x2;
184	cmat[count][3]=y2;
185	count++;
186	}
187	else if (fabs(s0-s2)<TOL&& /* Contour on side n2 -> n0 --- */
188	fabs(s2-cval)<TOL){
189	cmat[count][0]=x2;
190	cmat[count][1]=y2;
191	cmat[count][2]=x0;
192	cmat[count][3]=y0;
193	count++;
194	}
195	else if (fabs(s0-s2)<TOL){ /* Contour over entire element - */
196	cmat[count][0]=x0;
197	cmat[count][1]=y0;
198	cmat[count][2]=x1;
199	cmat[count][3]=y1;
200	count++;
201	cmat[count][0]=x1;
202	cmat[count][1]=y1;
203	cmat[count][2]=x2;
204	cmat[count][3]=y2;
205	count++;
206	cmat[count][0]=x2;
207	cmat[count][1]=y2;
208	cmat[count][2]=x0;
209	cmat[count][3]=y0;
210	count++;
211	}
212	else { /* Contour within element ----- */
213	fac=(cval-s0)/(s2-s0);
214	xa=x0+(x2-x0)*fac;
215	ya=y0+(y2-y0)*fac;
216	if(cval<s1){
217	if(s0!=s1) fac=(cval-s0)/(s1-s0);
218	else fac=1.0;
219	xb=x0+(x1-x0)*fac;
220	yb=y0+(y1-y0)*fac;
221	}
222	else{
223	if(s1!=s2) fac=(cval-s1)/(s2-s1);
224	else fac=1.0;
225	xb=x1+(x2-x1)*fac;
226	yb=y1+(y2-y1)*fac;
227	}
228	cmat[count][0]=xa;
229	cmat[count][1]=ya;
230	cmat[count][2]=xb;
231	cmat[count][3]=yb;
232	count++;
233	}
234	L10: continue;
235	}
236	*cnt=count;
237	return;
238	}
239
240
241
242
243
244

Note: See TracBrowser for help on using the browser.

root/gliderproc/trunk/MATLAB/opnml/MEX/contmex5.c

Download in other formats: