/[escript]/temp/finley/src/Assemble_jacobeans.c
ViewVC logotype

Annotation of /temp/finley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 776 - (hide annotations)
Wed Jul 12 00:07:31 2006 UTC (13 years, 3 months ago) by gross
Original Path: trunk/finley/src/Assemble_jacobeans.c
File MIME type: text/plain
File size: 17351 byte(s)
basic code for persistence of jacobeans added. routines are not called yet
1 jgs 150 /*
2 elspeth 616 ************************************************************
3     * Copyright 2006 by ACcESS MNRF *
4     * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7     * Licensed under the Open Software License version 3.0 *
8     * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11 jgs 150 */
12 jgs 82
13 jgs 150
14 jgs 82 /**************************************************************/
15    
16 jgs 150 /* Author: gross@access.edu.au */
17     /* Version: $Id$ */
18 jgs 82
19     /**************************************************************/
20    
21     #include "Assemble.h"
22     #include "Util.h"
23     #ifdef _OPENMP
24     #include <omp.h>
25     #endif
26    
27     /**************************************************************/
28 gross 776 /* */
29     /* Jacobean 1D */
30     /* */
31     void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,
32     dim_t numShape, dim_t numElements,index_t* nodes,
33     double* DSDv, dim_t numTest,double* DTDv,
34     double* dTdX, double* volume, index_t* element_id) {
35     #define DIM 1
36     #define LOCALDIM 1
37     register int e,q,s;
38     char error_msg[LenErrorMsg_MAX];
39     #pragma omp parallel
40     {
41     register double D,invD;
42     double X0[numShape];
43     #pragma omp for private(e,q,s) schedule(static)
44     for(e=0;e<numElements;e++){
45     for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
46     for (q=0;q<numQuad;q++) {
47     D=0;
48     for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
49     if (D==0.) {
50     sprintf(error_msg,"Assemble_jacobeans_1D: element %d (id %d) has length zero.",e,element_id[e]);
51     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
52     } else {
53     invD=1./D;
54     for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*invD;
55     volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
56     }
57     }
58     }
59 jgs 82
60 gross 776 }
61     #undef DIM
62     #undef LOCALDIM
63 jgs 82
64 gross 776 }
65     /**************************************************************/
66     /* */
67     /* Jacobean 2D */
68     /* */
69     void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
70     dim_t numShape, dim_t numElements,index_t* nodes,
71     double* DSDv, dim_t numTest,double* DTDv,
72     double* dTdX, double* volume, index_t* element_id) {
73     #define DIM 2
74     #define LOCALDIM 2
75     register int e,q,s;
76     char error_msg[LenErrorMsg_MAX];
77     #pragma omp parallel
78     {
79     register double dXdv00,dXdv10,dXdv01,dXdv11,
80     dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
81     double X0[numShape], X1[numShape];
82     #pragma omp for private(e,q,s) schedule(static)
83     for(e=0;e<numElements;e++){
84     for (s=0;s<numShape; s++) {
85     X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
86     X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
87     }
88     for (q=0;q<numQuad;q++) {
89     dXdv00=0;
90     dXdv10=0;
91     dXdv01=0;
92     dXdv11=0;
93     for (s=0;s<numShape; s++) {
94     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
95     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
96     dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
97     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
98     }
99     D = dXdv00*dXdv11 - dXdv01*dXdv10;
100     if (D==0.) {
101     sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
102     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
103     } else {
104     invD=1./D;
105     dvdX00= dXdv11*invD;
106     dvdX10=-dXdv10*invD;
107     dvdX01=-dXdv01*invD;
108     dvdX11= dXdv00*invD;
109 jgs 82
110 gross 776 for (s=0;s<numTest; s++) {
111     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
112     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
113     }
114     volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
115     }
116     }
117     }
118 jgs 82
119 gross 776 }
120     #undef DIM
121     #undef LOCALDIM
122     }
123     /**************************************************************/
124     /* */
125     /* Jacobean 3D */
126     /* */
127     void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
128     dim_t numShape, dim_t numElements,index_t* nodes,
129     double* DSDv, dim_t numTest,double* DTDv,
130     double* dTdX, double* volume, index_t* element_id) {
131     #define DIM 3
132     #define LOCALDIM 3
133     register int e,q,s;
134     char error_msg[LenErrorMsg_MAX];
135     #pragma omp parallel
136     {
137     register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
138     dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
139     double X0[numShape], X1[numShape], X2[numShape];
140     #pragma omp for private(e,q,s) schedule(static)
141     for(e=0;e<numElements;e++){
142     for (s=0;s<numShape; s++) {
143     X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
144     X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
145     X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
146     }
147     for (q=0;q<numQuad;q++) {
148     dXdv00=0;
149     dXdv10=0;
150     dXdv20=0;
151     dXdv01=0;
152     dXdv11=0;
153     dXdv21=0;
154     dXdv02=0;
155     dXdv12=0;
156     dXdv22=0;
157     for (s=0;s<numShape; s++) {
158     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
159     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
160     dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
161     dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
162     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
163     dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
164     dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
165     dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
166     dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
167     }
168     D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
169     if (D==0.) {
170     sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
171     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
172     } else {
173     invD=1./D;
174     dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
175     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
176     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
177     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
178     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
179     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
180     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
181     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
182     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
183 jgs 82
184 gross 776 for (s=0;s<numTest; s++) {
185     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX20;
186     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX21;
187     dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX22;
188     }
189     volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
190     }
191     }
192     }
193 jgs 82
194 gross 776 }
195     #undef DIM
196     #undef LOCALDIM
197     }
198     /**************************************************************/
199     /* */
200     /* Jacobean 2D manifold in 3D */
201     /* */
202     void Assemble_jacobeans_3D_M2D(double* coordinates, dim_t numQuad,double* QuadWeights,
203     dim_t numShape, dim_t numElements,index_t* nodes,
204     double* DSDv, dim_t numTest,double* DTDv,
205     double* dTdX, double* volume, index_t* element_id) {
206     #define DIM 3
207     #define LOCALDIM 2
208     register int e,q,s;
209     char error_msg[LenErrorMsg_MAX];
210     #pragma omp parallel
211     {
212     register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
213     dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;
214     double X0[numShape], X1[numShape], X2[numShape];
215     #pragma omp for private(e,q,s) schedule(static)
216     for(e=0;e<numElements;e++){
217     for (s=0;s<numShape; s++) {
218     X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
219     X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
220     X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
221     }
222     for (q=0;q<numQuad;q++) {
223     dXdv00=0;
224     dXdv10=0;
225     dXdv20=0;
226     dXdv01=0;
227     dXdv11=0;
228     dXdv21=0;
229     for (s=0;s<numShape; s++) {
230     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
231     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
232     dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
233     dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
234     dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
235     dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
236     }
237     m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
238     m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
239     m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
240     D=m00*m11-m01*m01;
241     if (D==0.) {
242     sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
243     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
244     } else {
245     invD=1./D;
246     dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
247     dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
248     dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
249     dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
250     dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
251     dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
252     for (s=0;s<numTest; s++) {
253     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
254     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
255     dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12;
256     }
257     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
258     }
259     }
260     }
261 jgs 82
262 gross 776 }
263     #undef DIM
264     #undef LOCALDIM
265     }
266     /**************************************************************/
267     /* */
268     /* Jacobean 1D manifold in 3D */
269     /* */
270     void Assemble_jacobeans_3D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,
271     dim_t numShape, dim_t numElements,index_t* nodes,
272     double* DSDv, dim_t numTest,double* DTDv,
273     double* dTdX, double* volume, index_t* element_id) {
274     #define DIM 3
275     #define LOCALDIM 1
276     register int e,q,s;
277     char error_msg[LenErrorMsg_MAX];
278     #pragma omp parallel
279     {
280     register double dXdv00,dXdv10,dXdv20,dvdX00,dvdX01,dvdX02,D,invD;
281     double X0[numShape], X1[numShape], X2[numShape];
282     #pragma omp for private(e,q,s) schedule(static)
283     for(e=0;e<numElements;e++){
284     for (s=0;s<numShape; s++) {
285     X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
286     X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
287     X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
288     }
289     for (q=0;q<numQuad;q++) {
290     dXdv00=0;
291     dXdv10=0;
292     dXdv20=0;
293     for (s=0;s<numShape; s++) {
294     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
295     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
296     dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
297     }
298     D=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
299     if (D==0.) {
300     sprintf(error_msg,"Assemble_jacobeans_3D_M1D: element %d (id %d) has length zero.",e,element_id[e]);
301     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
302     } else {
303     invD=1./D;
304     dvdX00=dXdv00*invD;
305     dvdX01=dXdv10*invD;
306     dvdX02=dXdv20*invD;
307     for (s=0;s<numTest; s++) {
308     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;
309     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;
310     dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02;
311     }
312     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
313     }
314     }
315 jgs 82 }
316    
317 gross 776 }
318     #undef DIM
319     #undef LOCALDIM
320     }
321     /**************************************************************/
322     /* */
323     /* Jacobean 1D manifold in 2D */
324     /* */
325     void Assemble_jacobeans_2D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,
326     dim_t numShape, dim_t numElements,index_t* nodes,
327     double* DSDv, dim_t numTest,double* DTDv,
328     double* dTdX, double* volume, index_t* element_id) {
329     #define DIM 2
330     #define LOCALDIM 1
331     register int e,q,s;
332     char error_msg[LenErrorMsg_MAX];
333     #pragma omp parallel
334 jgs 82 {
335 gross 776 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
336     double X0[numShape], X1[numShape];
337     #pragma omp for private(e,q,s) schedule(static)
338     for(e=0;e<numElements;e++){
339     for (s=0;s<numShape; s++) {
340     X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
341     X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
342     }
343     for (q=0;q<numQuad;q++) {
344     dXdv00=0;
345     dXdv10=0;
346     for (s=0;s<numShape; s++) {
347     dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
348     dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
349     }
350     D=dXdv00*dXdv00+dXdv10*dXdv10;
351     if (D==0.) {
352     sprintf(error_msg,"Assemble_jacobeans_2D_M1D: element %d (id %d) has length zero.",e,element_id[e]);
353     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
354     } else {
355     invD=1./D;
356     dvdX00=dXdv00*invD;
357     dvdX01=dXdv10*invD;
358     for (s=0;s<numTest; s++) {
359     dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;
360     dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;
361     }
362     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
363     }
364     }
365     }
366 jgs 82
367     }
368 gross 776 #undef DIM
369     #undef LOCALDIM
370 jgs 82 }
371     /*
372     * $Log$
373 jgs 150 *
374 jgs 82 */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26