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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 44248 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26