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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 781 - (show annotations)
Fri Jul 14 08:47:38 2006 UTC (13 years, 4 months ago) by gross
Original Path: trunk/finley/src/Assemble_jacobeans.c
File MIME type: text/plain
File size: 42353 byte(s)
grad functions linked into the persistent jacobean scheme
1 /*
2 ************************************************************
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 */
12
13
14 /**************************************************************/
15
16 /* Author: gross@access.edu.au */
17 /* Version: $Id$ */
18
19 /**************************************************************/
20
21 #include "Assemble.h"
22 #include "Util.h"
23 #ifdef _OPENMP
24 #include <omp.h>
25 #endif
26
27
28
29 /**************************************************************/
30 /* */
31 /* Jacobean 1D */
32 /* */
33 void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,
34 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
35 double* DSDv, dim_t numTest, double* DTDv,
36 double* dTdX, double* volume, index_t* element_id) {
37 #define DIM 1
38 #define LOCDIM 1
39 register int e,q,s;
40 char error_msg[LenErrorMsg_MAX];
41 #pragma omp parallel
42 {
43 register double D,invD;
44 double X0[numShape];
45 #pragma omp for private(e,q,s) schedule(static)
46 for(e=0;e<numElements;e++){
47 for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
48 for (q=0;q<numQuad;q++) {
49 D=0;
50 for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
51 if (D==0.) {
52 sprintf(error_msg,"Assemble_jacobeans_1D: element %d (id %d) has length zero.",e,element_id[e]);
53 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
54 } else {
55 invD=1./D;
56 for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*invD;
57 }
58 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
59 }
60 }
61
62 }
63 #undef DIM
64 #undef LOCDIM
65
66 }
67 /**************************************************************/
68 /* */
69 /* Jacobean 2D with area element */
70 /* */
71 void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
72 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
73 double* DSDv, dim_t numTest, double* DTDv,
74 double* dTdX, double* volume, index_t* element_id) {
75 #define DIM 2
76 #define LOCDIM 2
77 register int e,q,s;
78 char error_msg[LenErrorMsg_MAX];
79 #pragma omp parallel
80 {
81 register double dXdv00,dXdv10,dXdv01,dXdv11,
82 dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
83 double X0[numShape], X1[numShape];
84 #pragma omp for private(e,q,s) schedule(static)
85 for(e=0;e<numElements;e++){
86 for (s=0;s<numShape; s++) {
87 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
88 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
89 }
90 for (q=0;q<numQuad;q++) {
91 dXdv00=0;
92 dXdv10=0;
93 dXdv01=0;
94 dXdv11=0;
95 for (s=0;s<numShape; s++) {
96 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
97 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
98 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
99 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
100 }
101 D = dXdv00*dXdv11 - dXdv01*dXdv10;
102 if (D==0.) {
103 sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
104 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
105 } else {
106 invD=1./D;
107 dvdX00= dXdv11*invD;
108 dvdX10=-dXdv10*invD;
109 dvdX01=-dXdv01*invD;
110 dvdX11= dXdv00*invD;
111
112 for (s=0;s<numTest; s++) {
113 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10;
114 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11;
115 }
116 }
117 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
118 }
119 }
120
121 }
122 #undef DIM
123 #undef LOCDIM
124 }
125 /**************************************************************/
126 /* */
127 /* Jacobean 1D manifold in 2D and 1D elements */
128 /* */
129 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,double* QuadWeights,
130 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
131 double* DSDv, dim_t numTest, double* DTDv,
132 double* dTdX, double* volume, index_t* element_id) {
133 #define DIM 2
134 #define LOCDIM 1
135 register int e,q,s;
136 char error_msg[LenErrorMsg_MAX];
137 #pragma omp parallel
138 {
139 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
140 double X0[numShape], X1[numShape];
141 #pragma omp for private(e,q,s) schedule(static)
142 for(e=0;e<numElements;e++){
143 for (s=0;s<numShape; s++) {
144 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
145 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
146 }
147 for (q=0;q<numQuad;q++) {
148 dXdv00=0;
149 dXdv10=0;
150 for (s=0;s<numShape; s++) {
151 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
152 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
153 }
154 D=dXdv00*dXdv00+dXdv10*dXdv10;
155 if (D==0.) {
156 sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
157 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
158 } else {
159 invD=1./D;
160 dvdX00=dXdv00*invD;
161 dvdX01=dXdv10*invD;
162 for (s=0;s<numTest; s++) {
163 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00;
164 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01;
165 }
166 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
167 }
168 }
169 }
170
171 }
172 #undef DIM
173 #undef LOCDIM
174 }
175 /**************************************************************/
176 /* */
177 /* Jacobean 1D manifold in 2D and 1D elements woth contact */
178 /* */
179 void Assemble_jacobeans_2D_M1D_E1D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
180 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
181 double* DSDv, dim_t numTest, double* DTDv,
182 double* dTdX, double* volume, index_t* element_id) {
183 #define DIM 2
184 #define LOCDIM 1
185 register int e,q,s;
186 char error_msg[LenErrorMsg_MAX];
187 #pragma omp parallel
188 {
189 register double dXdv00_0,dXdv10_0,dvdX00_0,dvdX01_0,D_0,invD_0;
190 register double dXdv00_1,dXdv10_1,dvdX00_1,dvdX01_1,D_1,invD_1;
191 double X0[2*numShape], X1[2*numShape];
192 #pragma omp for private(e,q,s) schedule(static)
193 for(e=0;e<numElements;e++){
194 for (s=0;s<2*numShape; s++) {
195 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
196 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
197 }
198 for (q=0;q<numQuad;q++) {
199 dXdv00_0=0;
200 dXdv10_0=0;
201 dXdv00_1=0;
202 dXdv10_1=0;
203 for (s=0;s<numShape; s++) {
204 dXdv00_0+=X0[s] *DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
205 dXdv10_0+=X1[s] *DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
206 dXdv00_1+=X0[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
207 dXdv10_1+=X1[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
208 }
209 D_0=dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0;
210 D_1=dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1;
211 if (D_0 == 0.i || D_1 == 0.) {
212 sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
213 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
214 } else {
215 invD_0=1./D_0;
216 dvdX00_0=dXdv00_0*invD_0;
217 dvdX01_0=dXdv10_0*invD_0;
218 dvdX00_1=dXdv00_1*invD_1;
219 dvdX01_1=dXdv10_1*invD_1;
220 for (s=0;s<numTest; s++) {
221 dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0;
222 dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0;
223 dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1;
224 dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1;
225 }
226 volume[INDEX2(q,e,numQuad)]=(sqrt(D_0)+sqrt(D_1))/2.*QuadWeights[q];
227 }
228 }
229 }
230
231 }
232 #undef DIM
233 #undef LOCDIM
234 }
235 /**************************************************************/
236 /* */
237 /* Jacobean 1D manifold in 2D and 2D elements */
238 /* */
239 void Assemble_jacobeans_2D_M1D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
240 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
241 double* DSDv, dim_t numTest,double* DTDv,
242 double* dTdX, double* volume, index_t* element_id) {
243 #define DIM 2
244 #define LOCDIM 2
245 register int e,q,s;
246 char error_msg[LenErrorMsg_MAX];
247 #pragma omp parallel
248 {
249 register double dXdv00,dXdv10,dXdv01,dXdv11,
250 dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
251 double X0[numShape], X1[numShape];
252 #pragma omp for private(e,q,s) schedule(static)
253 for(e=0;e<numElements;e++){
254 for (s=0;s<numShape; s++) {
255 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
256 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
257 }
258 for (q=0;q<numQuad;q++) {
259 dXdv00=0;
260 dXdv10=0;
261 dXdv01=0;
262 dXdv11=0;
263 for (s=0;s<numShape; s++) {
264 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
265 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
266 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
267 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
268 }
269 D = dXdv00*dXdv11 - dXdv01*dXdv10;
270 if (D==0.) {
271 sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);
272 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
273 } else {
274 invD=1./D;
275 dvdX00= dXdv11*invD;
276 dvdX10=-dXdv10*invD;
277 dvdX01=-dXdv01*invD;
278 dvdX11= dXdv00*invD;
279
280 for (s=0;s<numTest; s++) {
281 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10;
282 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11;
283 }
284 }
285 volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];
286 }
287 }
288
289 }
290 #undef DIM
291 #undef LOCDIM
292 }
293 /**************************************************************/
294 /* */
295 /* Jacobean 1D manifold in 2D and 2D elements with contact */
296 /* */
297 void Assemble_jacobeans_2D_M1D_E2D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
298 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
299 double* DSDv, dim_t numTest,double* DTDv,
300 double* dTdX, double* volume, index_t* element_id) {
301 #define DIM 2
302 #define LOCDIM 2
303 register int e,q,s;
304 char error_msg[LenErrorMsg_MAX];
305 #pragma omp parallel
306 {
307 register double dXdv00_0,dXdv10_0,dXdv01_0,dXdv11_0,dvdX00_0,dvdX10_0,dvdX01_0,dvdX11_0, D_0,invD_0,
308 dXdv00_1,dXdv10_1,dXdv01_1,dXdv11_1,dvdX00_1,dvdX10_1,dvdX01_1,dvdX11_1, D_1,invD_1;
309 double X0[2*numShape], X1[2*numShape];
310 #pragma omp for private(e,q,s) schedule(static)
311 for(e=0;e<numElements;e++){
312 for (s=0;s<2*numShape; s++) {
313 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
314 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
315 }
316 for (q=0;q<numQuad;q++) {
317 dXdv00_0=0;
318 dXdv10_0=0;
319 dXdv01_0=0;
320 dXdv11_0=0;
321 dXdv00_1=0;
322 dXdv10_1=0;
323 dXdv01_1=0;
324 dXdv11_1=0;
325 for (s=0;s<numShape; s++) {
326 dXdv00_0+=X0[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
327 dXdv10_0+=X1[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
328 dXdv01_0+=X0[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
329 dXdv11_0+=X1[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
330 dXdv00_1+=X0[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
331 dXdv10_1+=X1[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
332 dXdv01_1+=X0[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
333 dXdv11_1+=X1[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
334 }
335 D_0 = dXdv00_0*dXdv11_0 - dXdv01_0*dXdv10_0;
336 D_1 = dXdv00_1*dXdv11_1 - dXdv01_1*dXdv10_1;
337 if ( (D_0 ==0.) || (D_1 ==0.) ) {
338 sprintf(error_msg,"Assemble_jacobeans_2D_E2D_C: element %d (id %d) has area zero.",e,element_id[e]);
339 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
340 } else {
341 invD_0=1./D_0;
342 dvdX00_0= dXdv11_0*invD_0;
343 dvdX10_0=-dXdv10_0*invD_0;
344 dvdX01_0=-dXdv01_0*invD_0;
345 dvdX11_0= dXdv00_0*invD_0;
346 invD_1=1./D_1;
347 dvdX00_1= dXdv11_1*invD_1;
348 dvdX10_1=-dXdv10_1*invD_1;
349 dvdX01_1=-dXdv01_1*invD_1;
350 dvdX11_1= dXdv00_1*invD_1;
351 for (s=0;s<numTest; s++) {
352 dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=
353 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0;
354 dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=
355 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0;
356 dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=
357 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1;
358 dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=
359 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1;
360 }
361 }
362 volume[INDEX2(q,e,numQuad)]=(sqrt(dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0)+sqrt(dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1))/2.*QuadWeights[q];
363 }
364 }
365
366 }
367 #undef DIM
368 #undef LOCDIM
369 }
370 /**************************************************************/
371 /* */
372 /* Jacobean 3D */
373 /* */
374 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
375 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
376 double* DSDv, dim_t numTest, double* DTDv,
377 double* dTdX, double* volume, index_t* element_id) {
378 #define DIM 3
379 #define LOCDIM 3
380 register int e,q,s;
381 char error_msg[LenErrorMsg_MAX];
382 #pragma omp parallel
383 {
384 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
385 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
386 double X0[numShape], X1[numShape], X2[numShape];
387 #pragma omp for private(e,q,s) schedule(static)
388 for(e=0;e<numElements;e++){
389 for (s=0;s<numShape; s++) {
390 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
391 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
392 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
393 }
394 for (q=0;q<numQuad;q++) {
395 dXdv00=0;
396 dXdv10=0;
397 dXdv20=0;
398 dXdv01=0;
399 dXdv11=0;
400 dXdv21=0;
401 dXdv02=0;
402 dXdv12=0;
403 dXdv22=0;
404 for (s=0;s<numShape; s++) {
405 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
406 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
407 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
408 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
409 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
410 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
411 dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
412 dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
413 dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
414 }
415 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
416 if (D==0.) {
417 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
418 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
419 } else {
420 invD=1./D;
421 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
422 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
423 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
424 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
425 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
426 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
427 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
428 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
429 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
430
431 for (s=0;s<numTest; s++) {
432 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
433 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;
434 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;
435 }
436 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
437 }
438 }
439 }
440
441 }
442 #undef DIM
443 #undef LOCDIM
444 }
445 /**************************************************************/
446 /* */
447 /* Jacobean 2D manifold in 3D with 3D elements */
448 /* */
449 void Assemble_jacobeans_3D_M2D_E3D(double* coordinates, dim_t numQuad,double* QuadWeights,
450 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
451 double* DSDv, dim_t numTest,double* DTDv,
452 double* dTdX, double* volume, index_t* element_id) {
453 #define DIM 3
454 #define LOCDIM 3
455 register int e,q,s;
456 char error_msg[LenErrorMsg_MAX];
457 #pragma omp parallel
458 {
459 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22, m0, m1, m2,
460 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
461 double X0[numShape], X1[numShape], X2[numShape];
462 #pragma omp for private(e,q,s) schedule(static)
463 for(e=0;e<numElements;e++){
464 for (s=0;s<numShape; s++) {
465 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
466 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
467 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
468 }
469 for (q=0;q<numQuad;q++) {
470 dXdv00=0;
471 dXdv10=0;
472 dXdv20=0;
473 dXdv01=0;
474 dXdv11=0;
475 dXdv21=0;
476 dXdv02=0;
477 dXdv12=0;
478 dXdv22=0;
479 for (s=0;s<numShape; s++) {
480 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
481 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
482 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
483 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
484 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
485 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
486 dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
487 dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
488 dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
489 }
490 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
491 if (D==0.) {
492 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
493 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
494 } else {
495 invD=1./D;
496 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
497 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
498 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
499 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
500 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
501 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
502 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
503 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
504 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
505
506 for (s=0;s<numTest; s++) {
507 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=
508 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
509 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=
510 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;
511 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=
512 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;
513 }
514 }
515 m0=dXdv10*dXdv21-dXdv20*dXdv11;
516 m1=dXdv20*dXdv01-dXdv00*dXdv21;
517 m2=dXdv00*dXdv11-dXdv10*dXdv01;
518 volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];
519 }
520 }
521
522 }
523 #undef DIM
524 #undef LOCDIM
525 }
526 /**************************************************************/
527 /* */
528 /* Jacobean 2D manifold in 3D with 3D elements on contact */
529 /* */
530 void Assemble_jacobeans_3D_M2D_E3D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
531 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
532 double* DSDv, dim_t numTest,double* DTDv,
533 double* dTdX, double* volume, index_t* element_id) {
534 #define DIM 3
535 #define LOCDIM 3
536 register int e,q,s;
537 char error_msg[LenErrorMsg_MAX];
538 #pragma omp parallel
539 {
540 register double dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,dXdv02_0,dXdv12_0,dXdv22_0, m0_0, m1_0, m2_0,
541 dvdX00_0,dvdX10_0,dvdX20_0,dvdX01_0,dvdX11_0,dvdX21_0,dvdX02_0,dvdX12_0,dvdX22_0, D_0,invD_0,
542 dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,dXdv02_1,dXdv12_1,dXdv22_1, m0_1, m1_1, m2_1,
543 dvdX00_1,dvdX10_1,dvdX20_1,dvdX01_1,dvdX11_1,dvdX21_1,dvdX02_1,dvdX12_1,dvdX22_1, D_1,invD_1;
544 double X0[2*numShape], X1[2*numShape], X2[2*numShape];
545 #pragma omp for private(e,q,s) schedule(static)
546 for(e=0;e<numElements;e++){
547 for (s=0;s<2*numShape; s++) {
548 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
549 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
550 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
551 }
552 for (q=0;q<numQuad;q++) {
553 dXdv00_0=0;
554 dXdv10_0=0;
555 dXdv20_0=0;
556 dXdv01_0=0;
557 dXdv11_0=0;
558 dXdv21_0=0;
559 dXdv02_0=0;
560 dXdv12_0=0;
561 dXdv22_0=0;
562 dXdv00_1=0;
563 dXdv10_1=0;
564 dXdv20_1=0;
565 dXdv01_1=0;
566 dXdv11_1=0;
567 dXdv21_1=0;
568 dXdv02_1=0;
569 dXdv12_1=0;
570 dXdv22_1=0;
571 for (s=0;s<numShape; s++) {
572 dXdv00_0+=X0[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
573 dXdv10_0+=X1[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
574 dXdv20_0+=X2[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
575 dXdv01_0+=X0[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
576 dXdv11_0+=X1[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
577 dXdv21_0+=X2[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
578 dXdv02_0+=X0[ s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
579 dXdv12_0+=X1[ s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
580 dXdv22_0+=X2[ s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
581 dXdv00_1+=X0[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
582 dXdv10_1+=X1[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
583 dXdv20_1+=X2[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
584 dXdv01_1+=X0[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
585 dXdv11_1+=X1[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
586 dXdv21_1+=X2[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
587 dXdv02_1+=X0[numShape+s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
588 dXdv12_1+=X1[numShape+s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
589 dXdv22_1+=X2[numShape+s]*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
590 }
591
592 D_0=dXdv00_0*(dXdv11_0*dXdv22_0-dXdv12_0*dXdv21_0)+dXdv01_0*(dXdv20_0*dXdv12_0-dXdv10_0*dXdv22_0)+dXdv02_0*(dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0);
593 D_1=dXdv00_1*(dXdv11_1*dXdv22_1-dXdv12_1*dXdv21_1)+dXdv01_1*(dXdv20_1*dXdv12_1-dXdv10_1*dXdv22_1)+dXdv02_1*(dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1);
594 if ( (D_0==0.) || (D_1 == 0.)) {
595 sprintf(error_msg,"Assemble_jacobeans_3D_C: element %d (id %d) has volume zero.",e,element_id[e]);
596 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
597 } else {
598 invD_0=1./D_0;
599 dvdX00_0=(dXdv11_0*dXdv22_0-dXdv12_0*dXdv21_0)*invD_0;
600 dvdX10_0=(dXdv20_0*dXdv12_0-dXdv10_0*dXdv22_0)*invD_0;
601 dvdX20_0=(dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0)*invD_0;
602 dvdX01_0=(dXdv02_0*dXdv21_0-dXdv01_0*dXdv22_0)*invD_0;
603 dvdX11_0=(dXdv00_0*dXdv22_0-dXdv20_0*dXdv02_0)*invD_0;
604 dvdX21_0=(dXdv01_0*dXdv20_0-dXdv00_0*dXdv21_0)*invD_0;
605 dvdX02_0=(dXdv01_0*dXdv12_0-dXdv02_0*dXdv11_0)*invD_0;
606 dvdX12_0=(dXdv02_0*dXdv10_0-dXdv00_0*dXdv12_0)*invD_0;
607 dvdX22_0=(dXdv00_0*dXdv11_0-dXdv01_0*dXdv10_0)*invD_0;
608 invD_1=1./D_1;
609 dvdX00_1=(dXdv11_1*dXdv22_1-dXdv12_1*dXdv21_1)*invD_1;
610 dvdX10_1=(dXdv20_1*dXdv12_1-dXdv10_1*dXdv22_1)*invD_1;
611 dvdX20_1=(dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1)*invD_1;
612 dvdX01_1=(dXdv02_1*dXdv21_1-dXdv01_1*dXdv22_1)*invD_1;
613 dvdX11_1=(dXdv00_1*dXdv22_1-dXdv20_1*dXdv02_1)*invD_1;
614 dvdX21_1=(dXdv01_1*dXdv20_1-dXdv00_1*dXdv21_1)*invD_1;
615 dvdX02_1=(dXdv01_1*dXdv12_1-dXdv02_1*dXdv11_1)*invD_1;
616 dvdX12_1=(dXdv02_1*dXdv10_1-dXdv00_1*dXdv12_1)*invD_1;
617 dvdX22_1=(dXdv00_1*dXdv11_1-dXdv01_1*dXdv10_1)*invD_1;
618
619 for (s=0;s<numTest; s++) {
620 dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=
621 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20_0;
622 dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=
623 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21_0;
624 dTdX[INDEX4( s,2,q,e,2*numTest,DIM,numQuad)]=
625 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22_0;
626 dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=
627 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20_1;
628 dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=
629 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21_1;
630 dTdX[INDEX4(numTest+s,2,q,e,2*numTest,DIM,numQuad)]=
631 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22_1;
632 }
633 }
634 m0_0=dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0;
635 m1_0=dXdv20_0*dXdv01_0-dXdv00_0*dXdv21_0;
636 m2_0=dXdv00_0*dXdv11_0-dXdv10_0*dXdv01_0;
637 m0_1=dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1;
638 m1_1=dXdv20_1*dXdv01_1-dXdv00_1*dXdv21_1;
639 m2_1=dXdv00_1*dXdv11_1-dXdv10_1*dXdv01_1;
640 volume[INDEX2(q,e,numQuad)]=(sqrt(m0_0*m0_0+m1_0*m1_0+m2_0*m2_0)+sqrt(m0_1*m0_1+m1_1*m1_1+m2_1*m2_1))/2.*QuadWeights[q];
641 }
642 }
643
644 }
645 #undef DIM
646 #undef LOCDIM
647 }
648 /**************************************************************/
649 /* */
650 /* Jacobean 2D manifold in 3D with 2D elements */
651 /* */
652 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
653 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
654 double* DSDv, dim_t numTest,double* DTDv,
655 double* dTdX, double* volume, index_t* element_id) {
656 #define DIM 3
657 #define LOCDIM 2
658 register int e,q,s;
659 char error_msg[LenErrorMsg_MAX];
660 #pragma omp parallel
661 {
662 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
663 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;
664 double X0[numShape], X1[numShape], X2[numShape];
665 #pragma omp for private(e,q,s) schedule(static)
666 for(e=0;e<numElements;e++){
667 for (s=0;s<numShape; s++) {
668 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
669 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
670 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
671 }
672 for (q=0;q<numQuad;q++) {
673 dXdv00=0;
674 dXdv10=0;
675 dXdv20=0;
676 dXdv01=0;
677 dXdv11=0;
678 dXdv21=0;
679 for (s=0;s<numShape; s++) {
680 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
681 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
682 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
683 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
684 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
685 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
686 }
687 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
688 m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
689 m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
690 D=m00*m11-m01*m01;
691 if (D==0.) {
692 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
693 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
694 } else {
695 invD=1./D;
696 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
697 dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
698 dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
699 dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
700 dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
701 dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
702 for (s=0;s<numTest; s++) {
703 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10;
704 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11;
705 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12;
706 }
707 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
708 }
709 }
710 }
711
712 }
713 #undef DIM
714 #undef LOCDIM
715 }
716 /**************************************************************/
717 /* */
718 /* Jacobean 2D manifold in 3D with 2D elements with contact */
719 /* */
720 void Assemble_jacobeans_3D_M2D_E2D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
721 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
722 double* DSDv, dim_t numTest,double* DTDv,
723 double* dTdX, double* volume, index_t* element_id) {
724 #define DIM 3
725 #define LOCDIM 2
726 register int e,q,s;
727 char error_msg[LenErrorMsg_MAX];
728 #pragma omp parallel
729 {
730 register double dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,m00_0,m01_0,m11_0,
731 dvdX00_0,dvdX01_0,dvdX02_0,dvdX10_0,dvdX11_0,dvdX12_0,D_0,invD_0,
732 dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,m00_1,m01_1,m11_1,
733 dvdX00_1,dvdX01_1,dvdX02_1,dvdX10_1,dvdX11_1,dvdX12_1,D_1,invD_1;
734 double X0[2*numShape], X1[2*numShape], X2[2*numShape];
735 #pragma omp for private(e,q,s) schedule(static)
736 for(e=0;e<numElements;e++){
737 for (s=0;s<2*numShape; s++) {
738 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
739 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
740 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
741 }
742 for (q=0;q<numQuad;q++) {
743 dXdv00_0=0;
744 dXdv10_0=0;
745 dXdv20_0=0;
746 dXdv01_0=0;
747 dXdv11_0=0;
748 dXdv21_0=0;
749 dXdv00_1=0;
750 dXdv10_1=0;
751 dXdv20_1=0;
752 dXdv01_1=0;
753 dXdv11_1=0;
754 dXdv21_1=0;
755 for (s=0;s<numShape; s++) {
756 dXdv00_0+=X0[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
757 dXdv10_0+=X1[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
758 dXdv20_0+=X2[ s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
759 dXdv01_0+=X0[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
760 dXdv11_0+=X1[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
761 dXdv21_0+=X2[ s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
762 dXdv00_1+=X0[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
763 dXdv10_1+=X1[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
764 dXdv20_1+=X2[numShape+s]*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
765 dXdv01_1+=X0[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
766 dXdv11_1+=X1[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
767 dXdv21_1+=X2[numShape+s]*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
768 }
769 m00_0=dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0+dXdv20_0*dXdv20_0;
770 m01_0=dXdv00_0*dXdv01_0+dXdv10_0*dXdv11_0+dXdv20_0*dXdv21_0;
771 m11_0=dXdv01_0*dXdv01_0+dXdv11_0*dXdv11_0+dXdv21_0*dXdv21_0;
772 D_0=m00_0*m11_0-m01_0*m01_0;
773 m00_1=dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1+dXdv20_1*dXdv20_1;
774 m01_1=dXdv00_1*dXdv01_1+dXdv10_1*dXdv11_1+dXdv20_1*dXdv21_1;
775 m11_1=dXdv01_1*dXdv01_1+dXdv11_1*dXdv11_1+dXdv21_1*dXdv21_1;
776 D_1=m00_1*m11_1-m01_1*m01_1;
777 if ( (D_0==0.) || (D_1 == 0.) ) {
778 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
779 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
780 } else {
781 invD_0=1./D_0;
782 dvdX00_0=( m00_0*dXdv00_0-m01_0*dXdv01_0)*invD_0;
783 dvdX01_0=( m00_0*dXdv10_0-m01_0*dXdv11_0)*invD_0;
784 dvdX02_0=( m00_0*dXdv20_0-m01_0*dXdv21_0)*invD_0;
785 dvdX10_0=(-m01_0*dXdv00_0+m11_0*dXdv01_0)*invD_0;
786 dvdX11_0=(-m01_0*dXdv10_0+m11_0*dXdv11_0)*invD_0;
787 dvdX12_0=(-m01_0*dXdv20_0+m11_0*dXdv21_0)*invD_0;
788 invD_1=1./D_1;
789 dvdX00_1=( m00_1*dXdv00_1-m01_1*dXdv01_1)*invD_1;
790 dvdX01_1=( m00_1*dXdv10_1-m01_1*dXdv11_1)*invD_1;
791 dvdX02_1=( m00_1*dXdv20_1-m01_1*dXdv21_1)*invD_1;
792 dvdX10_1=(-m01_1*dXdv00_1+m11_1*dXdv01_1)*invD_1;
793 dvdX11_1=(-m01_1*dXdv10_1+m11_1*dXdv11_1)*invD_1;
794 dvdX12_1=(-m01_1*dXdv20_1+m11_1*dXdv21_1)*invD_1;
795 for (s=0;s<numTest; s++) {
796 dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=
797 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0;
798 dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=
799 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0;
800 dTdX[INDEX4( s,2,q,e,2*numTest,DIM,numQuad)]=
801 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_0;
802 dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=
803 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1;
804 dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=
805 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1;
806 dTdX[INDEX4(numTest+s,2,q,e,2*numTest,DIM,numQuad)]=
807 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_1;
808 }
809 volume[INDEX2(q,e,numQuad)]=(sqrt(D_0)+sqrt(D_1))/2.*QuadWeights[q];
810 }
811 }
812 }
813
814 }
815 #undef DIM
816 #undef LOCDIM
817 }
818 /*
819 * $Log$
820 *
821 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26