NCCOOS Trac Projects: Top | Web | Platforms | Processing | Viz | Sprints | Sandbox | (Wind)

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

Revision 495 (checked in by cbc, 11 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 #define ELE(i,j,m) ele[i+m*j]
6
7
8
9 /* PROTOTYPES */
10 #ifdef __STDC__
11 void comp_bwidth(int,int *,int *);
12 void bandmsolve_c(int,
13                   double **,
14                   double *,
15                   int,
16                   int);
17 void grad_assem_lhs(int nn,
18           int ne,
19           int nbw,
20           double *x,
21           double *y,
22           int *ele,
23           double **lhs);
24 void grad_assem_rhs(int nn,
25           int ne,
26           int nbw,
27           double *x,
28           double *y,
29           int *ele,
30           double *q,
31           double *grax,
32           double *gray);
33
34 #endif
35
36 /************************************************************
37
38   ####     ##     #####  ######  #    #    ##     #   #
39  #    #   #  #      #    #       #    #   #  #     # #
40  #       #    #     #    #####   #    #  #    #     #
41  #  ###  ######     #    #       # ## #  ######     #
42  #    #  #    #     #    #       ##  ##  #    #     #
43   ####   #    #     #    ######  #    #  #    #     #
44
45 ************************************************************/
46
47 void mexFunction(int            nlhs,
48                  mxArray       *plhs[],
49                  int            nrhs,
50                  const mxArray *prhs[])
51 {
52    int i,j,nn,ne,nbw,nh,LHSout,NeedLHS;
53    double *x,*y,*q,*dele,*dqdx,*dqdy;
54    double **lhs;
55    int *ele;
56    double *temp,*LHS,*lhstemp;
57    
58 /* ---- gradmex5 will be called as :
59         [dqdx,dqdy,LHSout]=gradmex5(x,y,ele,q,LHS);
60                                       ------------------------------- */
61  /* ---- check I/O arguments ----------------------------------------- */
62    if (nrhs != 5)
63       mexErrMsgTxt("gradmex5 requires 5 input arguments.");
64    else if (nlhs != 3)
65       mexErrMsgTxt("gradmex5 requires 3 output arguments.");
66
67 /* ---- dereference input arrays ------------------------------------ */
68    x=mxGetPr(prhs[0]);
69    y=mxGetPr(prhs[1]);
70    dele=mxGetPr(prhs[2]);
71    q=mxGetPr(prhs[3]);
72    LHS=mxGetPr(prhs[4]);
73    nn=mxGetM(prhs[0]);
74    ne=mxGetM(prhs[2]);
75    fprintf(stderr,"NN=%d NE=%d\n",nn,ne);   
76
77    if (LHS[0]==1. & LHS[1]==0. & LHS[2]==1. & LHS[3]==1.)
78       NeedLHS=1;
79    else
80       NeedLHS=0;
81    
82 /* ---- allocate space for int representation of dele &
83         convert double element representation to int  &
84         shift node numbers toward 0 by 1 for proper indexing -------- */
85    ele=(int *)mxIvector(0,3*ne);
86    for (i=0;i<3*ne;i++){
87       ele[i]=(int)dele[i];
88       ele[i]=ele[i]-1;
89    }
90  
91 /* ---- get bandwidth ----------------------------------------------- */
92    comp_bwidth(ne,ele,&nbw);
93    nh=(nbw-1)/2; 
94    fprintf(stderr,"BW,HBW,NeedLHS = %d %d %d\n",nbw,nh,NeedLHS);
95
96 /* ---- allocate space for gradient lists  following
97         NRC allocation style                            ------------- */
98    dqdx = (double *)  mxDvector(0,nn-1);
99    dqdy = (double *)  mxDvector(0,nn-1);
100    
101 /* ---- Allocate, assemble, and triangularize LHS, if needed */   
102    lhs = (double **) mxDmatrix(0,nbw-1,0,nn-1);
103    if (NeedLHS){
104       grad_assem_lhs(nn,ne,nbw,x,y,ele,lhs);
105       bandmsolve_c(0,lhs,dqdx,nn,nh);
106       puts("LHS triangularized");
107    }
108    else{
109       for (i=0;i<nn;i++)
110          for (j=0;j<nbw;j++)
111             lhs[j][i]=LHS[i*nbw+j];       
112    }
113    
114    grad_assem_rhs(nn,ne,nbw,x,y,ele,q,dqdx,dqdy);
115    
116    puts("Backsolving for gradient");   
117    bandmsolve_c(1,lhs,dqdx,nn,nh);
118    bandmsolve_c(1,lhs,dqdy,nn,nh);
119    
120    puts("Done.");
121
122 /* ---- Set elements of return matricies,
123         pointed to by plhs[0] & plhs[1]
124                                         ----------------------------- */
125    plhs[0]=mxCreateDoubleMatrix(nn,1,mxREAL);
126    mxSetPr(plhs[0],dqdx);
127    plhs[1]=mxCreateDoubleMatrix(nn,1,mxREAL);
128    mxSetPr(plhs[1],dqdy);
129    plhs[2]=mxCreateDoubleMatrix(nbw,nn,mxREAL);
130    lhstemp=(double *)  mxDvector(0,nbw*nn-1);
131    for (i=0;i<nn;i++)
132          for (j=0;j<nbw;j++)
133             lhstemp[i*nbw+j]=lhs[j][i];
134            
135    mxSetPr(plhs[2],lhstemp);
136
137 /* --- No need to free memory allocated with "mxCalloc"; MATLAB
138    does this automatically.  The CMEX allocation functions in
139    "opnml_mex5_allocs.c" use mxCalloc. ------------------------------ */   
140    
141    return;
142
143 }
144
145
146 /*----------------------------------------------------------------------
147
148   ####   #####     ##    #####
149  #    #  #    #   #  #   #    #
150  #       #    #  #    #  #    #
151  #  ###  #####   ######  #    #
152  #    #  #   #   #    #  #    #
153   ####   #    #  #    #  #####
154  
155 ----------------------------------------------------------------------*/
156
157 void grad_assem_rhs(int nn,int ne,int nbw,
158           double *x,
159           double *y,
160           int *ele,
161           double *q,
162           double *grax,
163           double *gray)
164
165 {
166    double *ar,**av,**dx,**dy;
167    double area6,area12;
168    int i,j,l,m,n;
169
170 /* ---- ALLOCATE space for dx, dy, ar, av and dv -------------------- */
171    ar = (double *)  mxDvector(0,ne-1);
172    dx = (double **) mxDmatrix(0,2,0,ne-1);
173    dy = (double **) mxDmatrix(0,2,0,ne-1);   
174            
175 /* ---- Compute dx,dy and element areas ----------------------------- */
176    for(l=0;l<ne;l++){
177       dx[0][l]=x[ELE(l,1,ne)]-x[ELE(l,2,ne)];
178       dx[1][l]=x[ELE(l,2,ne)]-x[ELE(l,0,ne)];
179       dx[2][l]=x[ELE(l,0,ne)]-x[ELE(l,1,ne)];
180       dy[0][l]=y[ELE(l,1,ne)]-y[ELE(l,2,ne)];
181       dy[1][l]=y[ELE(l,2,ne)]-y[ELE(l,0,ne)];
182       dy[2][l]=y[ELE(l,0,ne)]-y[ELE(l,1,ne)];
183       ar[l]=.5*(x[ELE(l,0,ne)]*dy[0][l]+
184                 x[ELE(l,1,ne)]*dy[1][l]+
185                 x[ELE(l,2,ne)]*dy[2][l]);     
186    }
187    
188 /* ---- Assemble RHS (grax,gray) ------------------------------------ */
189    for(l=0;l<ne;l++){
190       for(j=0;j<3;j++){
191          m=ELE(l,j,ne);
192          for(i=0;i<3;i++){
193             grax[m]+=q[ELE(l,i,ne)]*dy[i][l]/(double)6.;
194             gray[m]-=q[ELE(l,i,ne)]*dx[i][l]/(double)6.;
195          }
196       }
197    }
198    puts("RHS (dqdx,dqdy) assembled");
199
200 }
201
202
203 void grad_assem_lhs(int nn,int ne,int nbw,
204           double *x,
205           double *y,
206           int *ele,
207           double **av)
208
209 {
210    double *ar,**dx,**dy;
211    double area6,area12;
212    int i,j,l,m,n,nh;
213
214    nh=(nbw-1)/2;
215    
216 /* ---- ALLOCATE local space for dx, dy, ar  -------------------- */
217    ar = (double *)  mxDvector(0,ne-1);
218    dx = (double **) mxDmatrix(0,2,0,ne-1);
219    dy = (double **) mxDmatrix(0,2,0,ne-1);   
220            
221 /* ---- Compute dx,dy and element areas ----------------------------- */
222    for(l=0;l<ne;l++){
223       dx[0][l]=x[ELE(l,1,ne)]-x[ELE(l,2,ne)];
224       dx[1][l]=x[ELE(l,2,ne)]-x[ELE(l,0,ne)];
225       dx[2][l]=x[ELE(l,0,ne)]-x[ELE(l,1,ne)];
226       dy[0][l]=y[ELE(l,1,ne)]-y[ELE(l,2,ne)];
227       dy[1][l]=y[ELE(l,2,ne)]-y[ELE(l,0,ne)];
228       dy[2][l]=y[ELE(l,0,ne)]-y[ELE(l,1,ne)];
229       ar[l]=.5*(x[ELE(l,0,ne)]*dy[0][l]+
230                 x[ELE(l,1,ne)]*dy[1][l]+
231                 x[ELE(l,2,ne)]*dy[2][l]);     
232    }
233    
234 /* ---- Assemble LHS (av) ------------------------------------------- */
235    for(l=0;l<ne;l++){
236       area6=ar[l]/(double)6.;
237       area12=ar[l]/(double)12.;
238       for(j=0;j<3;j++){
239          m=ELE(l,j,ne);
240          for(i=0;i<3;i++){
241             n=nh+ELE(l,i,ne)-ELE(l,j,ne);
242             if(i==j) av[n][m]+=area6;
243             else     av[n][m]+=area12;     
244          }
245       }
246    }
247    puts("LHS (av) assembled");
248    return;     
249 }
250
251 void comp_bwidth(int ne,int *ele,int *nbw)
252 {
253    int i, l01,l02,l12,nhm,nh;
254 /* ---- Compute half bandwidth -------------------------------------- */
255    for (i=0;i<ne;i++){
256       l01=abs(ELE(i,0,ne)-ELE(i,1,ne));
257       l02=abs(ELE(i,0,ne)-ELE(i,2,ne));
258       l12=abs(ELE(i,1,ne)-ELE(i,2,ne));
259       nhm=IMAX(l01,IMAX(l02,l12)); 
260       nh=IMAX(nh,nhm);   
261    }
262    *nbw=2*nh+1;
263 }
Note: See TracBrowser for help on using the browser.