/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_jacobians.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Assemble_jacobians.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 11 months ago) by caltinay
File size: 13738 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "Assemble.h"
18 #include "ShapeTable.h"
19 #include "Util.h"
20
21 // Unless the loops in here get complicated again, this file should be
22 // compiled with loop unrolling
23
24 /* input:
25
26 double* coordinates[DIM*(*)]
27 dim_t numQuad
28 double* QuadWeights[numQuad]
29 dim_t numShape
30 dim_t numElements
31 dim_t numNodes
32 index_t* nodes[numNodes*numElements] where NUMSIDES*numShape<=numNodes
33 double* DSDv[numShape*DIM*numQuad]
34 dim_t numTest
35 double* DTDv[LOCDIM*numTest*numQuad]
36 index_t* elementId[numElements]
37
38 output:
39
40 double* dTdX[DIM*numTest*NUMSIDES*numQuad*numElements]
41 double* volume[numQuad*numElements]
42
43 */
44
45 #define SCALING(_nsub_,_dim_) pow(1./(double)(_nsub_),1./(double)(_dim_))
46
47 namespace dudley {
48
49 /****************************************************************************/
50 //
51 // Jacobian 2D with area element
52 //
53 void Assemble_jacobians_2D(const double* coordinates, int numQuad,
54 dim_t numElements, int numNodes, const index_t* nodes,
55 double* dTdX, double* absD, double* quadweight,
56 const index_t* elementId)
57 {
58 const int DIM = 2;
59 const int numTest = 3; // hoping this is used in constant folding
60 *quadweight = (numQuad == 1) ? 1. / 2 : 1. / 6; // numQuad is 1 or 3
61 #pragma omp parallel for
62 for (index_t e = 0; e < numElements; e++) {
63 #define COMPDXDV0(P) coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
64 coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
65 coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
66
67 #define COMPDXDV1(P) coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
68 coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
69 coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
70
71 double dXdv00 = COMPDXDV0(0);
72 double dXdv10 = COMPDXDV0(1);
73 double dXdv01 = COMPDXDV1(0);
74 double dXdv11 = COMPDXDV1(1);
75 const double D = dXdv00 * dXdv11 - dXdv01 * dXdv10;
76 absD[e] = std::abs(D);
77 if (D == 0.) {
78 std::stringstream ss;
79 ss << "Assemble_jacobians_2D: element " << e
80 << " (id " << elementId[e] << ") has area zero.";
81 throw DudleyException(ss.str());
82 } else {
83 const double invD = 1. / D;
84 const double dvdX00 = dXdv11 * invD;
85 const double dvdX10 = -dXdv10 * invD;
86 const double dvdX01 = -dXdv01 * invD;
87 const double dvdX11 = dXdv00 * invD;
88 if (numQuad == 1) {
89 dTdX[INDEX4(0, 0, 0, e, numTest, DIM, numQuad)] =
90 DTDV_2D[0][0] * dvdX00 + DTDV_2D[1][1] * dvdX10;
91 dTdX[INDEX4(1, 0, 0, e, numTest, DIM, numQuad)] =
92 DTDV_2D[0][1] * dvdX00 + DTDV_2D[1][0] * dvdX10;
93 dTdX[INDEX4(2, 0, 0, e, numTest, DIM, numQuad)] =
94 DTDV_2D[2][0] * dvdX00 + DTDV_2D[2][1] * dvdX10;
95
96 dTdX[INDEX4(0, 1, 0, e, numTest, DIM, numQuad)] =
97 DTDV_2D[0][0] * dvdX01 + DTDV_2D[1][1] * dvdX11;
98 dTdX[INDEX4(1, 1, 0, e, numTest, DIM, numQuad)] =
99 DTDV_2D[0][1] * dvdX01 + DTDV_2D[1][0] * dvdX11;
100 dTdX[INDEX4(2, 1, 0, e, numTest, DIM, numQuad)] =
101 DTDV_2D[2][0] * dvdX01 + DTDV_2D[2][1] * dvdX11;
102
103 } else { // numQuad == 3
104 // relying on unroll loops to optimise this
105 for (int q = 0; q < numTest; ++q) {
106 dTdX[INDEX4(0, 0, q, e, numTest, DIM, numQuad)] =
107 DTDV_2D[0][0] * dvdX00 + DTDV_2D[1][1] * dvdX10;
108 dTdX[INDEX4(1, 0, q, e, numTest, DIM, numQuad)] =
109 DTDV_2D[0][1] * dvdX00 + DTDV_2D[1][0] * dvdX10;
110 dTdX[INDEX4(2, 0, q, e, numTest, DIM, numQuad)] =
111 DTDV_2D[2][0] * dvdX00 + DTDV_2D[2][1] * dvdX10;
112
113 dTdX[INDEX4(0, 1, q, e, numTest, DIM, numQuad)] =
114 DTDV_2D[0][0] * dvdX01 + DTDV_2D[1][1] * dvdX11;
115 dTdX[INDEX4(1, 1, q, e, numTest, DIM, numQuad)] =
116 DTDV_2D[0][1] * dvdX01 + DTDV_2D[1][0] * dvdX11;
117 dTdX[INDEX4(2, 1, q, e, numTest, DIM, numQuad)] =
118 DTDV_2D[2][0] * dvdX01 + DTDV_2D[2][1] * dvdX11;
119
120 }
121 }
122 }
123 } // end parallel for
124 #undef COMPDXDV0
125 #undef COMPDXDV1
126 }
127
128 //
129 // Jacobian 1D manifold in 2D and 1D elements
130 //
131 void Assemble_jacobians_2D_M1D_E1D(const double* coordinates, int numQuad,
132 dim_t numElements, int numNodes,
133 const index_t* nodes, double* dTdX,
134 double* absD, double* quadweight,
135 const index_t* elementId)
136 {
137 const int DIM = 2;
138 const int numTest = 2;
139 *quadweight = (numQuad == 1) ? 1.0 : 0.5; // numQuad is 1 or 2
140 #pragma omp parallel for
141 for (index_t e = 0; e < numElements; e++) {
142 double dXdv00 =
143 coordinates[INDEX2(0, nodes[INDEX2(0, e, numNodes)], DIM)] * (-1.) +
144 coordinates[INDEX2(0, nodes[INDEX2(1, e, numNodes)], DIM)];
145 double dXdv10 =
146 coordinates[INDEX2(1, nodes[INDEX2(0, e, numNodes)], DIM)] * (-1.) +
147 coordinates[INDEX2(1, nodes[INDEX2(1, e, numNodes)], DIM)];
148 const double D = dXdv00 * dXdv00 + dXdv10 * dXdv10;
149 if (D == 0.) {
150 std::stringstream ss;
151 ss << "Assemble_jacobians_2D_M1D_E1D: element " << e
152 << " (id " << elementId[e] << ") has length zero.";
153 throw DudleyException(ss.str());
154 } else {
155 const double invD = 1. / D;
156 const double dvdX00 = dXdv00 * invD;
157 const double dvdX01 = dXdv10 * invD;
158 // The number of quad points is 1 or 2
159 dTdX[INDEX4(0, 0, 0, e, numTest, DIM, numQuad)] = -1 * dvdX00;
160 dTdX[INDEX4(0, 1, 0, e, numTest, DIM, numQuad)] = -1 * dvdX01;
161 dTdX[INDEX4(1, 0, 0, e, numTest, DIM, numQuad)] = -1 * dvdX00;
162 dTdX[INDEX4(1, 1, 0, e, numTest, DIM, numQuad)] = -1 * dvdX01;
163 absD[e] = sqrt(D);
164 if (numQuad == 2) {
165 dTdX[INDEX4(0, 0, 1, e, numTest, DIM, numQuad)] = dvdX00;
166 dTdX[INDEX4(0, 1, 1, e, numTest, DIM, numQuad)] = dvdX01;
167 dTdX[INDEX4(1, 0, 1, e, numTest, DIM, numQuad)] = dvdX00;
168 dTdX[INDEX4(1, 1, 1, e, numTest, DIM, numQuad)] = dvdX01;
169 }
170 }
171 } // end parallel for
172 }
173
174 //
175 // Jacobian 3D
176 //
177 void Assemble_jacobians_3D(const double* coordinates, int numQuad,
178 dim_t numElements, int numNodes,
179 const index_t* nodes, double* dTdX, double* absD,
180 double* quadweight, const index_t* elementId)
181 {
182 const int DIM = 3;
183 const int numShape = 4;
184 const int numTest = 4;
185 *quadweight = (numQuad == 1) ? 1. / 6 : 1. / 24; // numQuad is 1 or 4
186
187 #pragma omp parallel for
188 for (index_t e = 0; e < numElements; e++) {
189 double dXdv00 = 0;
190 double dXdv10 = 0;
191 double dXdv20 = 0;
192 double dXdv01 = 0;
193 double dXdv11 = 0;
194 double dXdv21 = 0;
195 double dXdv02 = 0;
196 double dXdv12 = 0;
197 double dXdv22 = 0;
198 for (int s = 0; s < numShape; s++) {
199 const double X0_loc = coordinates[INDEX2(0, nodes[INDEX2(s, e, numNodes)], DIM)];
200 const double X1_loc = coordinates[INDEX2(1, nodes[INDEX2(s, e, numNodes)], DIM)];
201 const double X2_loc = coordinates[INDEX2(2, nodes[INDEX2(s, e, numNodes)], DIM)];
202 dXdv00 += X0_loc * DTDV_3D[s][0];
203 dXdv10 += X1_loc * DTDV_3D[s][0];
204 dXdv20 += X2_loc * DTDV_3D[s][0];
205 dXdv01 += X0_loc * DTDV_3D[s][1];
206 dXdv11 += X1_loc * DTDV_3D[s][1];
207 dXdv21 += X2_loc * DTDV_3D[s][1];
208 dXdv02 += X0_loc * DTDV_3D[s][2];
209 dXdv12 += X1_loc * DTDV_3D[s][2];
210 dXdv22 += X2_loc * DTDV_3D[s][2];
211 }
212 const double D = dXdv00 * (dXdv11 * dXdv22 - dXdv12 * dXdv21)
213 + dXdv01 * (dXdv20 * dXdv12 - dXdv10 * dXdv22)
214 + dXdv02 * (dXdv10 * dXdv21 - dXdv20 * dXdv11);
215 absD[e] = std::abs(D);
216 if (D == 0.) {
217 std::stringstream ss;
218 ss << "Assemble_jacobians_3D: element " << e
219 << " (id " << elementId[e] << ") has volume zero.";
220 throw DudleyException(ss.str());
221 } else {
222 const double invD = 1. / D;
223 const double dvdX00 = (dXdv11 * dXdv22 - dXdv12 * dXdv21) * invD;
224 const double dvdX10 = (dXdv20 * dXdv12 - dXdv10 * dXdv22) * invD;
225 const double dvdX20 = (dXdv10 * dXdv21 - dXdv20 * dXdv11) * invD;
226 const double dvdX01 = (dXdv02 * dXdv21 - dXdv01 * dXdv22) * invD;
227 const double dvdX11 = (dXdv00 * dXdv22 - dXdv20 * dXdv02) * invD;
228 const double dvdX21 = (dXdv01 * dXdv20 - dXdv00 * dXdv21) * invD;
229 const double dvdX02 = (dXdv01 * dXdv12 - dXdv02 * dXdv11) * invD;
230 const double dvdX12 = (dXdv02 * dXdv10 - dXdv00 * dXdv12) * invD;
231 const double dvdX22 = (dXdv00 * dXdv11 - dXdv01 * dXdv10) * invD;
232 for (int q = 0; q < numQuad; q++) {
233 for (int s = 0; s < numTest; s++) {
234 dTdX[INDEX4(s, 0, q, e, numTest, DIM, numQuad)] =
235 DTDV_3D[s][0] * dvdX00 + DTDV_3D[s][1] * dvdX10
236 + DTDV_3D[s][2] * dvdX20;
237 dTdX[INDEX4(s, 1, q, e, numTest, DIM, numQuad)] =
238 DTDV_3D[s][0] * dvdX01 + DTDV_3D[s][1] * dvdX11
239 + DTDV_3D[s][2] * dvdX21;
240 dTdX[INDEX4(s, 2, q, e, numTest, DIM, numQuad)] =
241 DTDV_3D[s][0] * dvdX02 + DTDV_3D[s][1] * dvdX12
242 + DTDV_3D[s][2] * dvdX22;
243 }
244 }
245 }
246 } // end parallel for
247 }
248
249 //
250 // Jacobian 2D manifold in 3D with 2D elements
251 //
252 void Assemble_jacobians_3D_M2D_E2D(const double* coordinates, int numQuad,
253 dim_t numElements, int numNodes,
254 const index_t* nodes, double* dTdX,
255 double* absD, double* quadweight,
256 const index_t* elementId)
257 {
258 const int DIM = 3;
259 const double DTDV[3][2] = { {-1., -1.}, {1., 0.}, {0., 1.} };
260 const int numShape = 3;
261 const int numTest = 3;
262 *quadweight = (numQuad == 1) ? 1. / 2 : 1. / 6; // numQuad is 1 or 3
263 #pragma omp parallel for
264 for (index_t e = 0; e < numElements; e++) {
265 double dXdv00 = 0;
266 double dXdv10 = 0;
267 double dXdv20 = 0;
268 double dXdv01 = 0;
269 double dXdv11 = 0;
270 double dXdv21 = 0;
271 for (int s = 0; s < numShape; s++) {
272 const double X0_loc = coordinates[INDEX2(0, nodes[INDEX2(s, e, numNodes)], DIM)];
273 const double X1_loc = coordinates[INDEX2(1, nodes[INDEX2(s, e, numNodes)], DIM)];
274 const double X2_loc = coordinates[INDEX2(2, nodes[INDEX2(s, e, numNodes)], DIM)];
275 dXdv00 += X0_loc * DTDV[s][0];
276 dXdv10 += X1_loc * DTDV[s][0];
277 dXdv20 += X2_loc * DTDV[s][0];
278 dXdv01 += X0_loc * DTDV[s][1];
279 dXdv11 += X1_loc * DTDV[s][1];
280 dXdv21 += X2_loc * DTDV[s][1];
281 }
282 const double m00 = dXdv00 * dXdv00 + dXdv10 * dXdv10 + dXdv20 * dXdv20;
283 const double m01 = dXdv00 * dXdv01 + dXdv10 * dXdv11 + dXdv20 * dXdv21;
284 const double m11 = dXdv01 * dXdv01 + dXdv11 * dXdv11 + dXdv21 * dXdv21;
285 const double D = m00 * m11 - m01 * m01;
286 absD[e] = sqrt(D);
287 if (D == 0.) {
288 std::stringstream ss;
289 ss << "Assemble_jacobians_3D_M2D: element " << e
290 << " (id " << elementId[e] << ") has area zero.";
291 throw DudleyException(ss.str());
292 } else {
293 const double invD = 1. / D;
294 const double dvdX00 = (m00 * dXdv00 - m01 * dXdv01) * invD;
295 const double dvdX01 = (m00 * dXdv10 - m01 * dXdv11) * invD;
296 const double dvdX02 = (m00 * dXdv20 - m01 * dXdv21) * invD;
297 const double dvdX10 = (-m01 * dXdv00 + m11 * dXdv01) * invD;
298 const double dvdX11 = (-m01 * dXdv10 + m11 * dXdv11) * invD;
299 const double dvdX12 = (-m01 * dXdv20 + m11 * dXdv21) * invD;
300 for (int q = 0; q < numQuad; q++) {
301 for (int s = 0; s < numTest; s++) {
302 dTdX[INDEX4(s, 0, q, e, numTest, DIM, numQuad)] =
303 DTDV[s][0] * dvdX00 + DTDV[s][1] * dvdX10;
304 dTdX[INDEX4(s, 1, q, e, numTest, DIM, numQuad)] =
305 DTDV[s][0] * dvdX01 + DTDV[s][1] * dvdX11;
306 dTdX[INDEX4(s, 2, q, e, numTest, DIM, numQuad)] =
307 DTDV[s][0] * dvdX02 + DTDV[s][1] * dvdX12;
308 }
309 }
310 }
311 } // end parallel for
312 }
313
314 } // namespace dudley
315

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/dudley/src/Assemble_jacobeans.cpp:5567-5588 /branches/4.0fordebian/dudley/src/Assemble_jacobians.cpp:5567-5588 /branches/complex/dudley/src/Assemble_jacobeans.cpp:5866-5937 /branches/complex/dudley/src/Assemble_jacobians.cpp:5866-5937 /branches/diaplayground/dudley/src/Assemble_jacobians.cpp:4940-5147 /branches/lapack2681/finley/src/Assemble_jacobeans.cpp:2682-2741 /branches/lapack2681/finley/src/Assemble_jacobians.cpp:2682-2741 /branches/pasowrap/dudley/src/Assemble_jacobeans.cpp:3661-3674 /branches/pasowrap/dudley/src/Assemble_jacobians.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Assemble_jacobeans.cpp:3871-3891 /branches/py3_attempt2/dudley/src/Assemble_jacobians.cpp:3871-3891 /branches/restext/finley/src/Assemble_jacobeans.cpp:2610-2624 /branches/restext/finley/src/Assemble_jacobians.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Assemble_jacobeans.cpp:3669-3791 /branches/ripleygmg_from_3668/dudley/src/Assemble_jacobians.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_jacobeans.cpp:2569-2590 /branches/stage3.0/finley/src/Assemble_jacobians.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Assemble_jacobeans.cpp:3471-3974 /branches/symbolic_from_3470/dudley/src/Assemble_jacobians.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Assemble_jacobeans.cpp:3517-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Assemble_jacobians.cpp:3517-3974 /release/3.0/finley/src/Assemble_jacobeans.cpp:2591-2601 /release/3.0/finley/src/Assemble_jacobians.cpp:2591-2601 /release/4.0/dudley/src/Assemble_jacobeans.cpp:5380-5406 /release/4.0/dudley/src/Assemble_jacobians.cpp:5380-5406 /trunk/dudley/src/Assemble_jacobeans.cpp:4257-4344 /trunk/dudley/src/Assemble_jacobians.cpp:5898-5962,5982-6007 /trunk/ripley/test/python/dudley/src/Assemble_jacobeans.cpp:3480-3515 /trunk/ripley/test/python/dudley/src/Assemble_jacobians.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26