/[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 1628 - (show annotations)
Fri Jul 11 13:12:46 2008 UTC (11 years, 5 months ago) by phornby
File MIME type: text/plain
File size: 44283 byte(s)

Merge in /branches/windows_from_1456_trunk_1620_merged_in branch.

You will find a preserved pre-merge trunk in tags under tags/trunk_at_1625.
That will be useful for diffing & checking on my stupidity.

Here is a list of the conflicts and their resolution at this
point in time.


=================================================================================
(LLWS == looks like white space).

finley/src/Assemble_addToSystemMatrix.c - resolve to branch - unused var. may be wrong.....
finley/src/CPPAdapter/SystemMatrixAdapter.cpp - resolve to branch - LLWS
finley/src/CPPAdapter/MeshAdapter.cpp - resolve to branch - LLWS
paso/src/PCG.c - resolve to branch - unused var fixes.
paso/src/SolverFCT.c - resolve to branch - LLWS
paso/src/FGMRES.c - resolve to branch - LLWS
paso/src/Common.h - resolve to trunk version. It's omp.h's include... not sure it's needed,
but for the sake of saftey.....
paso/src/Functions.c - resolve to branch version, indentation/tab removal and return error
on bad unimplemented Paso_FunctionCall.
paso/src/SolverFCT_solve.c - resolve to branch version, unused vars
paso/src/SparseMatrix_MatrixVector.c - resolve to branch version, unused vars.
escript/src/Utils.cpp - resloved to branch, needs WinSock2.h
escript/src/DataExpanded.cpp - resolved to branch version - LLWS
escript/src/DataFactory.cpp - resolve to branch version
=================================================================================

This currently passes tests on linux (debian), but is not checked on windows or Altix yet.

This checkin is to make a trunk I can check out for windows to do tests on it.

Known outstanding problem is in the operator=() method of exceptions
causing warning messages on the intel compilers.

May the God of doughnuts have mercy on my soul.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26