/[escript]/branches/domexper/dudley/src/Assemble_jacobeans.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26