OOFEM 3.0
Loading...
Searching...
No Matches
fei3dhexatriquad.C
Go to the documentation of this file.
1/*
2 *
3 * ##### ##### ###### ###### ### ###
4 * ## ## ## ## ## ## ## ### ##
5 * ## ## ## ## #### #### ## # ##
6 * ## ## ## ## ## ## ## ##
7 * ## ## ## ## ## ## ## ##
8 * ##### ##### ## ###### ## ##
9 *
10 *
11 * OOFEM : Object Oriented Finite Element Code
12 *
13 * Copyright (C) 1993 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "fei3dhexatriquad.h"
36#include "intarray.h"
37#include "floatarray.h"
38#include "floatarrayf.h"
39#include "floatmatrix.h"
40#include "floatmatrixf.h"
41#include "gaussintegrationrule.h"
42
43namespace oofem {
44
45FloatArrayF<27>
46FEI3dHexaTriQuad :: evalN(const FloatArrayF<3> &lcoords)
47{
48 //auto [u, v, w] = lcoords;
49 double u = lcoords[0];
50 double v = lcoords[1];
51 double w = lcoords[2];
52
53 std::array<double, 3> a = {0.5 * ( u - 1.0 ) * u, 0.5 * ( u + 1.0 ) * u, 1.0 - u * u};
54 std::array<double, 3> b = {0.5 * ( v - 1.0 ) * v, 0.5 * ( v + 1.0 ) * v, 1.0 - v * v};
55 std::array<double, 3> c = {0.5 * ( w - 1.0 ) * w, 0.5 * ( w + 1.0 ) * w, 1.0 - w * w};
56
57 return {
58 a [ 0 ] * b [ 0 ] * c [ 1 ],
59 a [ 0 ] * b [ 1 ] * c [ 1 ],
60 a [ 1 ] * b [ 1 ] * c [ 1 ],
61 a [ 1 ] * b [ 0 ] * c [ 1 ],
62 a [ 0 ] * b [ 0 ] * c [ 0 ],
63 a [ 0 ] * b [ 1 ] * c [ 0 ],
64 a [ 1 ] * b [ 1 ] * c [ 0 ],
65 a [ 1 ] * b [ 0 ] * c [ 0 ],
66 a [ 0 ] * b [ 2 ] * c [ 1 ],
67 a [ 2 ] * b [ 1 ] * c [ 1 ],
68 a [ 1 ] * b [ 2 ] * c [ 1 ],
69 a [ 2 ] * b [ 0 ] * c [ 1 ],
70 a [ 0 ] * b [ 2 ] * c [ 0 ],
71 a [ 2 ] * b [ 1 ] * c [ 0 ],
72 a [ 1 ] * b [ 2 ] * c [ 0 ],
73 a [ 2 ] * b [ 0 ] * c [ 0 ],
74 a [ 0 ] * b [ 0 ] * c [ 2 ],
75 a [ 0 ] * b [ 1 ] * c [ 2 ],
76 a [ 1 ] * b [ 1 ] * c [ 2 ],
77 a [ 1 ] * b [ 0 ] * c [ 2 ],
78 a [ 2 ] * b [ 2 ] * c [ 1 ],
79 a [ 0 ] * b [ 2 ] * c [ 2 ],
80 a [ 2 ] * b [ 2 ] * c [ 0 ],
81 a [ 2 ] * b [ 1 ] * c [ 2 ],
82 a [ 1 ] * b [ 2 ] * c [ 2 ],
83 a [ 2 ] * b [ 0 ] * c [ 2 ],
84 a [ 2 ] * b [ 2 ] * c [ 2 ]
85 };
86}
87
88void
89FEI3dHexaTriQuad :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
90{
91#if 0
92 answer = evalN(lcoords);
93#else
94 double u = lcoords.at(1);
95 double v = lcoords.at(2);
96 double w = lcoords.at(3);
97
98 double a[] = {
99 0.5 * ( u - 1.0 ) * u, 0.5 * ( u + 1.0 ) * u, 1.0 - u * u
100 };
101 double b[] = {
102 0.5 * ( v - 1.0 ) * v, 0.5 * ( v + 1.0 ) * v, 1.0 - v * v
103 };
104 double c[] = {
105 0.5 * ( w - 1.0 ) * w, 0.5 * ( w + 1.0 ) * w, 1.0 - w * w
106 };
107
108 answer.resize(27);
109
110 answer.at(5) = a [ 0 ] * b [ 0 ] * c [ 0 ];
111 answer.at(8) = a [ 1 ] * b [ 0 ] * c [ 0 ];
112 answer.at(16) = a [ 2 ] * b [ 0 ] * c [ 0 ];
113
114 answer.at(6) = a [ 0 ] * b [ 1 ] * c [ 0 ];
115 answer.at(7) = a [ 1 ] * b [ 1 ] * c [ 0 ];
116 answer.at(14) = a [ 2 ] * b [ 1 ] * c [ 0 ];
117
118 answer.at(13) = a [ 0 ] * b [ 2 ] * c [ 0 ];
119 answer.at(15) = a [ 1 ] * b [ 2 ] * c [ 0 ];
120 answer.at(22) = a [ 2 ] * b [ 2 ] * c [ 0 ];
121
122 answer.at(1) = a [ 0 ] * b [ 0 ] * c [ 1 ];
123 answer.at(4) = a [ 1 ] * b [ 0 ] * c [ 1 ];
124 answer.at(12) = a [ 2 ] * b [ 0 ] * c [ 1 ];
125
126 answer.at(2) = a [ 0 ] * b [ 1 ] * c [ 1 ];
127 answer.at(3) = a [ 1 ] * b [ 1 ] * c [ 1 ];
128 answer.at(10) = a [ 2 ] * b [ 1 ] * c [ 1 ];
129
130 answer.at(9) = a [ 0 ] * b [ 2 ] * c [ 1 ];
131 answer.at(11) = a [ 1 ] * b [ 2 ] * c [ 1 ];
132 answer.at(21) = a [ 2 ] * b [ 2 ] * c [ 1 ];
133
134 answer.at(17) = a [ 0 ] * b [ 0 ] * c [ 2 ];
135 answer.at(20) = a [ 1 ] * b [ 0 ] * c [ 2 ];
136 answer.at(26) = a [ 2 ] * b [ 0 ] * c [ 2 ];
137
138 answer.at(18) = a [ 0 ] * b [ 1 ] * c [ 2 ];
139 answer.at(19) = a [ 1 ] * b [ 1 ] * c [ 2 ];
140 answer.at(24) = a [ 2 ] * b [ 1 ] * c [ 2 ];
141
142 answer.at(23) = a [ 0 ] * b [ 2 ] * c [ 2 ];
143 answer.at(25) = a [ 1 ] * b [ 2 ] * c [ 2 ];
144 answer.at(27) = a [ 2 ] * b [ 2 ] * c [ 2 ];
145#endif
146}
147
148
149void
150FEI3dHexaTriQuad :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
151{
152 double u = lcoords.at(1);
153 double v = lcoords.at(2);
154
155 double a[] = {
156 0.5 * ( u - 1. ) * u, 0.5 * ( u + 1. ) * u, 1. - u * u
157 };
158 double b[] = {
159 0.5 * ( v - 1. ) * v, 0.5 * ( v + 1. ) * v, 1. - v * v
160 };
161
162 answer.resize(9);
163
164 answer.at(1) = a [ 0 ] * b [ 0 ];
165 answer.at(2) = a [ 1 ] * b [ 0 ];
166 answer.at(5) = a [ 2 ] * b [ 0 ];
167
168 answer.at(4) = a [ 0 ] * b [ 1 ];
169 answer.at(3) = a [ 1 ] * b [ 1 ];
170 answer.at(7) = a [ 2 ] * b [ 1 ];
171
172 answer.at(8) = a [ 0 ] * b [ 2 ];
173 answer.at(6) = a [ 1 ] * b [ 2 ];
174 answer.at(9) = a [ 2 ] * b [ 2 ];
175}
176
177
178double
179FEI3dHexaTriQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
180{
181 FloatArray e1, e2, dNdu(9), dNdv(9);
182
183 double u = lcoords.at(1);
184 double v = lcoords.at(2);
185
186 double a[] = {
187 0.5 * ( u - 1. ) * u, 0.5 * ( u + 1. ) * u, 1. - u * u
188 };
189 double b[] = {
190 0.5 * ( v - 1. ) * v, 0.5 * ( v + 1. ) * v, 1. - v * v
191 };
192
193 double da[] = {
194 u - 0.5, u + 0.5, -2. * u
195 };
196 double db[] = {
197 v - 0.5, v + 0.5, -2. * v
198 };
199
200 dNdu.at(1) = da [ 0 ] * b [ 0 ];
201 dNdu.at(2) = da [ 1 ] * b [ 0 ];
202 dNdu.at(5) = da [ 2 ] * b [ 0 ];
203
204 dNdu.at(4) = da [ 0 ] * b [ 1 ];
205 dNdu.at(3) = da [ 1 ] * b [ 1 ];
206 dNdu.at(7) = da [ 2 ] * b [ 1 ];
207
208 dNdu.at(8) = da [ 0 ] * b [ 2 ];
209 dNdu.at(6) = da [ 1 ] * b [ 2 ];
210 dNdu.at(9) = da [ 2 ] * b [ 2 ];
211
212
213 dNdv.at(1) = a [ 0 ] * db [ 0 ];
214 dNdv.at(2) = a [ 1 ] * db [ 0 ];
215 dNdv.at(5) = a [ 2 ] * db [ 0 ];
216
217 dNdv.at(4) = a [ 0 ] * db [ 1 ];
218 dNdv.at(3) = a [ 1 ] * db [ 1 ];
219 dNdv.at(7) = a [ 2 ] * db [ 1 ];
220
221 dNdv.at(8) = a [ 0 ] * db [ 2 ];
222 dNdv.at(6) = a [ 1 ] * db [ 2 ];
223 dNdv.at(9) = a [ 2 ] * db [ 2 ];
224
225 const auto &snodes = this->computeLocalSurfaceMapping(isurf);
226 for ( int i = 1; i <= 9; ++i ) {
227 e1.add( dNdu.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
228 e2.add( dNdv.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
229 }
230
231 answer.beVectorProductOf(e1, e2);
232 return answer.normalize_giveNorm();
233}
234
235
236IntArray
237FEI3dHexaTriQuad :: computeLocalSurfaceMapping(int isurf) const
238{
239 if ( isurf == 1 ) {
240 return { 2, 1, 4, 3, 9, 12, 11, 10, 21};
241 } else if ( isurf == 2 ) {
242 return { 5, 6, 7, 8, 13, 14, 15, 16, 22};
243 } else if ( isurf == 3 ) {
244 return { 1, 2, 6, 5, 9, 18, 13, 17, 23};
245 } else if ( isurf == 4 ) {
246 return { 2, 3, 7, 6, 10, 19, 14, 18, 24};
247 } else if ( isurf == 5 ) {
248 return { 3, 4, 8, 7, 11, 20, 15, 19, 25};
249 } else if ( isurf == 6 ) {
250 return { 4, 1, 5, 8, 12, 17, 16, 20, 26};
251 } else {
252 throw std::range_error("invalid surface number");
253 }
254}
255
256
257FloatMatrixF<3,27>
258FEI3dHexaTriQuad :: evaldNdxi(const FloatArrayF<3> &lcoords)
259{
260 //auto [u, v, w] = lcoords;
261 double u = lcoords[0];
262 double v = lcoords[1];
263 double w = lcoords[2];
264
265 std::array<double, 3> a = {0.5 * ( u - 1.0 ) * u, 0.5 * ( u + 1.0 ) * u, 1.0 - u * u};
266 std::array<double, 3> b = {0.5 * ( v - 1.0 ) * v, 0.5 * ( v + 1.0 ) * v, 1.0 - v * v};
267 std::array<double, 3> c = {0.5 * ( w - 1.0 ) * w, 0.5 * ( w + 1.0 ) * w, 1.0 - w * w};
268
269 std::array<double, 3> da = {-0.5 + u, 0.5 + u, -2.0 * u};
270 std::array<double, 3> db = {-0.5 + v, 0.5 + v, -2.0 * v};
271 std::array<double, 3> dc = {-0.5 + w, 0.5 + w, -2.0 * w};
272
273 return {
274 da [ 0 ] * b [ 0 ] * c [ 1 ], a [ 0 ] * db [ 0 ] * c [ 1 ], a [ 0 ] * b [ 0 ] * dc [ 1 ],
275 da [ 0 ] * b [ 1 ] * c [ 1 ], a [ 0 ] * db [ 1 ] * c [ 1 ], a [ 0 ] * b [ 1 ] * dc [ 1 ],
276 da [ 1 ] * b [ 1 ] * c [ 1 ], a [ 1 ] * db [ 1 ] * c [ 1 ], a [ 1 ] * b [ 1 ] * dc [ 1 ],
277 da [ 1 ] * b [ 0 ] * c [ 1 ], a [ 1 ] * db [ 0 ] * c [ 1 ], a [ 1 ] * b [ 0 ] * dc [ 1 ],
278 da [ 0 ] * b [ 0 ] * c [ 0 ], a [ 0 ] * db [ 0 ] * c [ 0 ], a [ 0 ] * b [ 0 ] * dc [ 0 ],
279 da [ 0 ] * b [ 1 ] * c [ 0 ], a [ 0 ] * db [ 1 ] * c [ 0 ], a [ 0 ] * b [ 1 ] * dc [ 0 ],
280 da [ 1 ] * b [ 1 ] * c [ 0 ], a [ 1 ] * db [ 1 ] * c [ 0 ], a [ 1 ] * b [ 1 ] * dc [ 0 ],
281 da [ 1 ] * b [ 0 ] * c [ 0 ], a [ 1 ] * db [ 0 ] * c [ 0 ], a [ 1 ] * b [ 0 ] * dc [ 0 ],
282 da [ 0 ] * b [ 2 ] * c [ 1 ], a [ 0 ] * db [ 2 ] * c [ 1 ], a [ 0 ] * b [ 2 ] * dc [ 1 ],
283 da [ 2 ] * b [ 1 ] * c [ 1 ], a [ 2 ] * db [ 1 ] * c [ 1 ], a [ 2 ] * b [ 1 ] * dc [ 1 ],
284 da [ 1 ] * b [ 2 ] * c [ 1 ], a [ 1 ] * db [ 2 ] * c [ 1 ], a [ 1 ] * b [ 2 ] * dc [ 1 ],
285 da [ 2 ] * b [ 0 ] * c [ 1 ], a [ 2 ] * db [ 0 ] * c [ 1 ], a [ 2 ] * b [ 0 ] * dc [ 1 ],
286 da [ 0 ] * b [ 2 ] * c [ 0 ], a [ 0 ] * db [ 2 ] * c [ 0 ], a [ 0 ] * b [ 2 ] * dc [ 0 ],
287 da [ 2 ] * b [ 1 ] * c [ 0 ], a [ 2 ] * db [ 1 ] * c [ 0 ], a [ 2 ] * b [ 1 ] * dc [ 0 ],
288 da [ 1 ] * b [ 2 ] * c [ 0 ], a [ 1 ] * db [ 2 ] * c [ 0 ], a [ 1 ] * b [ 2 ] * dc [ 0 ],
289 da [ 2 ] * b [ 0 ] * c [ 0 ], a [ 2 ] * db [ 0 ] * c [ 0 ], a [ 2 ] * b [ 0 ] * dc [ 0 ],
290 da [ 0 ] * b [ 0 ] * c [ 2 ], a [ 0 ] * db [ 0 ] * c [ 2 ], a [ 0 ] * b [ 0 ] * dc [ 2 ],
291 da [ 0 ] * b [ 1 ] * c [ 2 ], a [ 0 ] * db [ 1 ] * c [ 2 ], a [ 0 ] * b [ 1 ] * dc [ 2 ],
292 da [ 1 ] * b [ 1 ] * c [ 2 ], a [ 1 ] * db [ 1 ] * c [ 2 ], a [ 1 ] * b [ 1 ] * dc [ 2 ],
293 da [ 1 ] * b [ 0 ] * c [ 2 ], a [ 1 ] * db [ 0 ] * c [ 2 ], a [ 1 ] * b [ 0 ] * dc [ 2 ],
294 da [ 2 ] * b [ 2 ] * c [ 1 ], a [ 2 ] * db [ 2 ] * c [ 1 ], a [ 2 ] * b [ 2 ] * dc [ 1 ],
295 da [ 2 ] * b [ 2 ] * c [ 0 ], a [ 2 ] * db [ 2 ] * c [ 0 ], a [ 2 ] * b [ 2 ] * dc [ 0 ],
296 da [ 0 ] * b [ 2 ] * c [ 2 ], a [ 0 ] * db [ 2 ] * c [ 2 ], a [ 0 ] * b [ 2 ] * dc [ 2 ],
297 da [ 2 ] * b [ 1 ] * c [ 2 ], a [ 2 ] * db [ 1 ] * c [ 2 ], a [ 2 ] * b [ 1 ] * dc [ 2 ],
298 da [ 1 ] * b [ 2 ] * c [ 2 ], a [ 1 ] * db [ 2 ] * c [ 2 ], a [ 1 ] * b [ 2 ] * dc [ 2 ],
299 da [ 2 ] * b [ 0 ] * c [ 2 ], a [ 2 ] * db [ 0 ] * c [ 2 ], a [ 2 ] * b [ 0 ] * dc [ 2 ],
300 da [ 2 ] * b [ 2 ] * c [ 2 ], a [ 2 ] * db [ 2 ] * c [ 2 ], a [ 2 ] * b [ 2 ] * dc [ 2 ],
301 };
302}
303
304
305void
306FEI3dHexaTriQuad :: evaldNdxi(FloatMatrix &dN, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
307{
308#if 0
309 dN = evaldNdxi(lcoords);
310#else
311 double u, v, w;
312 u = lcoords.at(1);
313 v = lcoords.at(2);
314 w = lcoords.at(3);
315
316 // Helpers expressions;
317 double a[] = {
318 0.5 * ( u - 1.0 ) * u, 0.5 * ( u + 1.0 ) * u, 1.0 - u * u
319 };
320 double b[] = {
321 0.5 * ( v - 1.0 ) * v, 0.5 * ( v + 1.0 ) * v, 1.0 - v * v
322 };
323 double c[] = {
324 0.5 * ( w - 1.0 ) * w, 0.5 * ( w + 1.0 ) * w, 1.0 - w * w
325 };
326
327 double da[] = {
328 -0.5 + u, 0.5 + u, -2.0 * u
329 };
330 double db[] = {
331 -0.5 + v, 0.5 + v, -2.0 * v
332 };
333 double dc[] = {
334 -0.5 + w, 0.5 + w, -2.0 * w
335 };
336
337 dN.resize(27, 3);
338
339 dN.at(5, 1) = da [ 0 ] * b [ 0 ] * c [ 0 ];
340 dN.at(8, 1) = da [ 1 ] * b [ 0 ] * c [ 0 ];
341 dN.at(16, 1) = da [ 2 ] * b [ 0 ] * c [ 0 ];
342
343 dN.at(6, 1) = da [ 0 ] * b [ 1 ] * c [ 0 ];
344 dN.at(7, 1) = da [ 1 ] * b [ 1 ] * c [ 0 ];
345 dN.at(14, 1) = da [ 2 ] * b [ 1 ] * c [ 0 ];
346
347 dN.at(13, 1) = da [ 0 ] * b [ 2 ] * c [ 0 ];
348 dN.at(15, 1) = da [ 1 ] * b [ 2 ] * c [ 0 ];
349 dN.at(22, 1) = da [ 2 ] * b [ 2 ] * c [ 0 ];
350
351 dN.at(1, 1) = da [ 0 ] * b [ 0 ] * c [ 1 ];
352 dN.at(4, 1) = da [ 1 ] * b [ 0 ] * c [ 1 ];
353 dN.at(12, 1) = da [ 2 ] * b [ 0 ] * c [ 1 ];
354
355 dN.at(2, 1) = da [ 0 ] * b [ 1 ] * c [ 1 ];
356 dN.at(3, 1) = da [ 1 ] * b [ 1 ] * c [ 1 ];
357 dN.at(10, 1) = da [ 2 ] * b [ 1 ] * c [ 1 ];
358
359 dN.at(9, 1) = da [ 0 ] * b [ 2 ] * c [ 1 ];
360 dN.at(11, 1) = da [ 1 ] * b [ 2 ] * c [ 1 ];
361 dN.at(21, 1) = da [ 2 ] * b [ 2 ] * c [ 1 ];
362
363 dN.at(17, 1) = da [ 0 ] * b [ 0 ] * c [ 2 ];
364 dN.at(20, 1) = da [ 1 ] * b [ 0 ] * c [ 2 ];
365 dN.at(26, 1) = da [ 2 ] * b [ 0 ] * c [ 2 ];
366
367 dN.at(18, 1) = da [ 0 ] * b [ 1 ] * c [ 2 ];
368 dN.at(19, 1) = da [ 1 ] * b [ 1 ] * c [ 2 ];
369 dN.at(24, 1) = da [ 2 ] * b [ 1 ] * c [ 2 ];
370
371 dN.at(23, 1) = da [ 0 ] * b [ 2 ] * c [ 2 ];
372 dN.at(25, 1) = da [ 1 ] * b [ 2 ] * c [ 2 ];
373 dN.at(27, 1) = da [ 2 ] * b [ 2 ] * c [ 2 ];
374
375 //
376
377 dN.at(5, 2) = a [ 0 ] * db [ 0 ] * c [ 0 ];
378 dN.at(8, 2) = a [ 1 ] * db [ 0 ] * c [ 0 ];
379 dN.at(16, 2) = a [ 2 ] * db [ 0 ] * c [ 0 ];
380
381 dN.at(6, 2) = a [ 0 ] * db [ 1 ] * c [ 0 ];
382 dN.at(7, 2) = a [ 1 ] * db [ 1 ] * c [ 0 ];
383 dN.at(14, 2) = a [ 2 ] * db [ 1 ] * c [ 0 ];
384
385 dN.at(13, 2) = a [ 0 ] * db [ 2 ] * c [ 0 ];
386 dN.at(15, 2) = a [ 1 ] * db [ 2 ] * c [ 0 ];
387 dN.at(22, 2) = a [ 2 ] * db [ 2 ] * c [ 0 ];
388
389 dN.at(1, 2) = a [ 0 ] * db [ 0 ] * c [ 1 ];
390 dN.at(4, 2) = a [ 1 ] * db [ 0 ] * c [ 1 ];
391 dN.at(12, 2) = a [ 2 ] * db [ 0 ] * c [ 1 ];
392
393 dN.at(2, 2) = a [ 0 ] * db [ 1 ] * c [ 1 ];
394 dN.at(3, 2) = a [ 1 ] * db [ 1 ] * c [ 1 ];
395 dN.at(10, 2) = a [ 2 ] * db [ 1 ] * c [ 1 ];
396
397 dN.at(9, 2) = a [ 0 ] * db [ 2 ] * c [ 1 ];
398 dN.at(11, 2) = a [ 1 ] * db [ 2 ] * c [ 1 ];
399 dN.at(21, 2) = a [ 2 ] * db [ 2 ] * c [ 1 ];
400
401 dN.at(17, 2) = a [ 0 ] * db [ 0 ] * c [ 2 ];
402 dN.at(20, 2) = a [ 1 ] * db [ 0 ] * c [ 2 ];
403 dN.at(26, 2) = a [ 2 ] * db [ 0 ] * c [ 2 ];
404
405 dN.at(18, 2) = a [ 0 ] * db [ 1 ] * c [ 2 ];
406 dN.at(19, 2) = a [ 1 ] * db [ 1 ] * c [ 2 ];
407 dN.at(24, 2) = a [ 2 ] * db [ 1 ] * c [ 2 ];
408
409 dN.at(23, 2) = a [ 0 ] * db [ 2 ] * c [ 2 ];
410 dN.at(25, 2) = a [ 1 ] * db [ 2 ] * c [ 2 ];
411 dN.at(27, 2) = a [ 2 ] * db [ 2 ] * c [ 2 ];
412
413 //
414
415 dN.at(5, 3) = a [ 0 ] * b [ 0 ] * dc [ 0 ];
416 dN.at(8, 3) = a [ 1 ] * b [ 0 ] * dc [ 0 ];
417 dN.at(16, 3) = a [ 2 ] * b [ 0 ] * dc [ 0 ];
418
419 dN.at(6, 3) = a [ 0 ] * b [ 1 ] * dc [ 0 ];
420 dN.at(7, 3) = a [ 1 ] * b [ 1 ] * dc [ 0 ];
421 dN.at(14, 3) = a [ 2 ] * b [ 1 ] * dc [ 0 ];
422
423 dN.at(13, 3) = a [ 0 ] * b [ 2 ] * dc [ 0 ];
424 dN.at(15, 3) = a [ 1 ] * b [ 2 ] * dc [ 0 ];
425 dN.at(22, 3) = a [ 2 ] * b [ 2 ] * dc [ 0 ];
426
427 dN.at(1, 3) = a [ 0 ] * b [ 0 ] * dc [ 1 ];
428 dN.at(4, 3) = a [ 1 ] * b [ 0 ] * dc [ 1 ];
429 dN.at(12, 3) = a [ 2 ] * b [ 0 ] * dc [ 1 ];
430
431 dN.at(2, 3) = a [ 0 ] * b [ 1 ] * dc [ 1 ];
432 dN.at(3, 3) = a [ 1 ] * b [ 1 ] * dc [ 1 ];
433 dN.at(10, 3) = a [ 2 ] * b [ 1 ] * dc [ 1 ];
434
435 dN.at(9, 3) = a [ 0 ] * b [ 2 ] * dc [ 1 ];
436 dN.at(11, 3) = a [ 1 ] * b [ 2 ] * dc [ 1 ];
437 dN.at(21, 3) = a [ 2 ] * b [ 2 ] * dc [ 1 ];
438
439 dN.at(17, 3) = a [ 0 ] * b [ 0 ] * dc [ 2 ];
440 dN.at(20, 3) = a [ 1 ] * b [ 0 ] * dc [ 2 ];
441 dN.at(26, 3) = a [ 2 ] * b [ 0 ] * dc [ 2 ];
442
443 dN.at(18, 3) = a [ 0 ] * b [ 1 ] * dc [ 2 ];
444 dN.at(19, 3) = a [ 1 ] * b [ 1 ] * dc [ 2 ];
445 dN.at(24, 3) = a [ 2 ] * b [ 1 ] * dc [ 2 ];
446
447 dN.at(23, 3) = a [ 0 ] * b [ 2 ] * dc [ 2 ];
448 dN.at(25, 3) = a [ 1 ] * b [ 2 ] * dc [ 2 ];
449 dN.at(27, 3) = a [ 2 ] * b [ 2 ] * dc [ 2 ];
450#endif
451}
452
453
454std::pair<double, FloatMatrixF<3,27>>
455FEI3dHexaTriQuad :: evaldNdx(const FloatArrayF<3> &lcoords, const FEICellGeometry &cellgeo)
456{
457 auto dNduvw = evaldNdxi(lcoords);
458 FloatMatrixF<3,27> coords;
459 for ( int i = 0; i < 27; i++ ) {
461 coords.setColumn(cellgeo.giveVertexCoordinates(i+1), i);
462 }
463 auto jacT = dotT(dNduvw, coords);
464 return {det(jacT), dot(inv(jacT), dNduvw)};
465}
466
467double
468FEI3dHexaTriQuad :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
469{
470 auto tmp = evaldNdx(lcoords, cellgeo);
471 answer = transpose(tmp.second);
472 return tmp.first;
473}
474
475double
476FEI3dHexaTriQuad :: evalNXIntegral(int iSurf, const FEICellGeometry &cellgeo) const
477{
478 const auto &fNodes = this->computeLocalSurfaceMapping(iSurf);
479
480 const auto &c1 = cellgeo.giveVertexCoordinates( fNodes.at(1) );
481 const auto &c2 = cellgeo.giveVertexCoordinates( fNodes.at(2) );
482 const auto &c3 = cellgeo.giveVertexCoordinates( fNodes.at(3) );
483 const auto &c4 = cellgeo.giveVertexCoordinates( fNodes.at(4) );
484 const auto &c5 = cellgeo.giveVertexCoordinates( fNodes.at(5) );
485 const auto &c6 = cellgeo.giveVertexCoordinates( fNodes.at(6) );
486 const auto &c7 = cellgeo.giveVertexCoordinates( fNodes.at(7) );
487 const auto &c8 = cellgeo.giveVertexCoordinates( fNodes.at(8) );
488 const auto &c9 = cellgeo.giveVertexCoordinates( fNodes.at(9) );
489
490 // Generated with Mathematica (rather unwieldy expression, tried to simplify it as good as possible, but it could probably be better)
491 return (
492 c1(2) * (
493 c2(1) * ( -3 * c3(0) - 3 * c4(0) + 18 * c6(0) - 4 * c7(0) + 18 * c8(0) + 24 * c9(0) ) +
494 c3(1) * ( 3 * c2(0) - 3 * c4(0) + 2 * c5(0) + 2 * c6(0) - 2 * c7(0) - 2 * c8(0) ) +
495 c4(1) * ( 3 * c2(0) + 3 * c3(0) - 18 * c5(0) + 4 * c6(0) - 18 * c7(0) - 24 * c9(0) ) +
496 c5(1) * ( -2 * c3(0) + 18 * c4(0) + 12 * c6(0) + 24 * c7(0) - 108 * c8(0) - 144 * c9(0) ) +
497 c6(1) * ( -18 * c2(0) - 2 * c3(0) - 4 * c4(0) - 12 * c5(0) - 4 * c7(0) + 24 * c8(0) + 16 * c9(0) ) +
498 c7(1) * ( 4 * c2(0) + 2 * c3(0) + 18 * c4(0) - 24 * c5(0) + 4 * c6(0) + 12 * c8(0) - 16 * c9(0) ) +
499 c8(1) * ( -18 * c2(0) + 2 * c3(0) + 108 * c5(0) - 24 * c6(0) - 12 * c7(0) + 144 * c9(0) ) +
500 c9(1) * ( -24 * c2(0) + 24 * c4(0) + 144 * c5(0) - 16 * c6(0) + 16 * c7(0) - 144 * c8(0) )
501 ) +
502 c2(2) * (
503 c1(1) * ( 3 * c3(0) + 3 * c4(0) - 18 * c6(0) + 4 * c7(0) - 18 * c8(0) - 24 * c9(0) ) +
504 c3(1) * ( -3 * c1(0) - 3 * c4(0) + 18 * c5(0) + 18 * c7(0) - 4 * c8(0) + 24 * c9(0) ) +
505 c4(1) * ( -3 * c1(0) + 3 * c3(0) - 2 * c5(0) + 2 * c6(0) + 2 * c7(0) - 2 * c8(0) ) +
506 c5(1) * ( -18 * c3(0) + 2 * c4(0) + 108 * c6(0) - 24 * c7(0) - 12 * c8(0) + 144 * c9(0) ) +
507 c6(1) * ( 18 * c1(0) - 2 * c4(0) - 108 * c5(0) + 12 * c7(0) + 24 * c8(0) - 144 * c9(0) ) +
508 c7(1) * ( -4 * c1(0) - 18 * c3(0) - 2 * c4(0) + 24 * c5(0) - 12 * c6(0) - 4 * c8(0) + 16 * c9(0) ) +
509 c8(1) * ( 18 * c1(0) + 4 * c3(0) + 2 * c4(0) + 12 * c5(0) - 24 * c6(0) + 4 * c7(0) - 16 * c9(0) ) +
510 c9(1) * ( 24 * c1(0) - 24 * c3(0) - 144 * c5(0) + 144 * c6(0) - 16 * c7(0) + 16 * c8(0) )
511 ) +
512 c3(2) * (
513 c1(1) * ( -3 * c2(0) + 3 * c4(0) - 2 * c5(0) - 2 * c6(0) + 2 * c7(0) + 2 * c8(0) ) +
514 c2(1) * ( 3 * c1(0) + 3 * c4(0) - 18 * c5(0) - 18 * c7(0) + 4 * c8(0) - 24 * c9(0) ) +
515 c4(1) * ( -3 * c1(0) - 3 * c2(0) - 4 * c5(0) + 18 * c6(0) + 18 * c8(0) + 24 * c9(0) ) +
516 c5(1) * ( 2 * c1(0) + 18 * c2(0) + 4 * c4(0) + 12 * c6(0) - 24 * c7(0) + 4 * c8(0) - 16 * c9(0) ) +
517 c6(1) * ( 2 * c1(0) - 18 * c4(0) - 12 * c5(0) + 108 * c7(0) - 24 * c8(0) + 144 * c9(0) ) +
518 c7(1) * ( -2 * c1(0) + 18 * c2(0) + 24 * c5(0) - 108 * c6(0) + 12 * c8(0) - 144 * c9(0) ) +
519 c8(1) * ( -2 * c1(0) - 4 * c2(0) - 18 * c4(0) - 4 * c5(0) + 24 * c6(0) - 12 * c7(0) + 16 * c9(0) ) +
520 c9(1) * ( 24 * c2(0) - 24 * c4(0) + 16 * c5(0) - 144 * c6(0) + 144 * c7(0) - 16 * c8(0) )
521 ) +
522 c4(2) * (
523 c1(1) * ( -3 * c2(0) - 3 * c3(0) + 18 * c5(0) - 4 * c6(0) + 18 * c7(0) + 24 * c9(0) ) +
524 c2(1) * ( 3 * c1(0) - 3 * c3(0) + 2 * c5(0) - 2 * c6(0) - 2 * c7(0) + 2 * c8(0) ) +
525 c3(1) * ( 3 * c1(0) + 3 * c2(0) + 4 * c5(0) - 18 * c6(0) - 18 * c8(0) - 24 * c9(0) ) +
526 c5(1) * ( -18 * c1(0) - 2 * c2(0) - 4 * c3(0) - 4 * c6(0) + 24 * c7(0) - 12 * c8(0) + 16 * c9(0) ) +
527 c6(1) * ( 4 * c1(0) + 2 * c2(0) + 18 * c3(0) + 4 * c5(0) + 12 * c7(0) - 24 * c8(0) - 16 * c9(0) ) +
528 c7(1) * ( -18 * c1(0) + 2 * c2(0) - 24 * c5(0) - 12 * c6(0) + 108 * c8(0) + 144 * c9(0) ) +
529 c8(1) * ( -2 * c2(0) + 18 * c3(0) + 12 * c5(0) + 24 * c6(0) - 108 * c7(0) - 144 * c9(0) ) +
530 c9(1) * ( -24 * c1(0) + 24 * c3(0) - 16 * c5(0) + 16 * c6(0) - 144 * c7(0) + 144 * c8(0) )
531 ) +
532 c5(2) * (
533 c1(1) * ( 2 * c3(0) - 18 * c4(0) - 12 * c6(0) - 24 * c7(0) + 108 * c8(0) + 144 * c9(0) ) +
534 c2(1) * ( 18 * c3(0) - 2 * c4(0) - 108 * c6(0) + 24 * c7(0) + 12 * c8(0) - 144 * c9(0) ) +
535 c3(1) * ( -2 * c1(0) - 18 * c2(0) - 4 * c4(0) - 12 * c6(0) + 24 * c7(0) - 4 * c8(0) + 16 * c9(0) ) +
536 c4(1) * ( 18 * c1(0) + 2 * c2(0) + 4 * c3(0) + 4 * c6(0) - 24 * c7(0) + 12 * c8(0) - 16 * c9(0) ) +
537 c6(1) * ( 12 * c1(0) + 108 * c2(0) + 12 * c3(0) - 4 * c4(0) + 32 * c7(0) + 32 * c8(0) - 192 * c9(0) ) +
538 c7(1) * ( 24 * c1(0) - 24 * c2(0) - 24 * c3(0) + 24 * c4(0) - 32 * c6(0) + 32 * c8(0) ) +
539 c8(1) * ( -108 * c1(0) - 12 * c2(0) + 4 * c3(0) - 12 * c4(0) - 32 * c6(0) - 32 * c7(0) + 192 * c9(0) ) +
540 c9(1) * ( -144 * c1(0) + 144 * c2(0) - 16 * c3(0) + 16 * c4(0) + 192 * c6(0) - 192 * c8(0) )
541 ) +
542 c6(2) * (
543 c1(1) * ( 18 * c2(0) + 2 * c3(0) + 4 * c4(0) + 12 * c5(0) + 4 * c7(0) - 24 * c8(0) - 16 * c9(0) ) +
544 c2(1) * ( -18 * c1(0) + 2 * c4(0) + 108 * c5(0) - 12 * c7(0) - 24 * c8(0) + 144 * c9(0) ) +
545 c3(1) * ( -2 * c1(0) + 18 * c4(0) + 12 * c5(0) - 108 * c7(0) + 24 * c8(0) - 144 * c9(0) ) +
546 c4(1) * ( -4 * c1(0) - 2 * c2(0) - 18 * c3(0) - 4 * c5(0) - 12 * c7(0) + 24 * c8(0) + 16 * c9(0) ) +
547 c5(1) * ( -12 * c1(0) - 108 * c2(0) - 12 * c3(0) + 4 * c4(0) - 32 * c7(0) - 32 * c8(0) + 192 * c9(0) ) +
548 c8(1) * ( 24 * c1(0) + 24 * c2(0) - 24 * c3(0) - 24 * c4(0) + 32 * c5(0) - 32 * c7(0) ) +
549 c7(1) * ( -4 * c1(0) + 12 * c2(0) + 108 * c3(0) + 12 * c4(0) + 32 * c5(0) + 32 * c8(0) - 192 * c9(0) ) +
550 c9(1) * ( 16 * c1(0) - 144 * c2(0) + 144 * c3(0) - 16 * c4(0) - 192 * c5(0) + 192 * c7(0) )
551 ) +
552 c7(2) * (
553 c1(1) * ( -4 * c2(0) - 2 * c3(0) - 18 * c4(0) + 24 * c5(0) - 4 * c6(0) - 12 * c8(0) + 16 * c9(0) ) +
554 c2(1) * ( 4 * c1(0) + 18 * c3(0) + 2 * c4(0) - 24 * c5(0) + 12 * c6(0) + 4 * c8(0) - 16 * c9(0) ) +
555 c3(1) * ( 2 * c1(0) - 18 * c2(0) - 24 * c5(0) + 108 * c6(0) - 12 * c8(0) + 144 * c9(0) ) +
556 c4(1) * ( 18 * c1(0) - 2 * c2(0) + 24 * c5(0) + 12 * c6(0) - 108 * c8(0) - 144 * c9(0) ) +
557 c5(1) * ( -24 * c1(0) + 24 * c2(0) + 24 * c3(0) - 24 * c4(0) + 32 * c6(0) - 32 * c8(0) ) +
558 c6(1) * ( 4 * c1(0) - 12 * c2(0) - 108 * c3(0) - 12 * c4(0) - 32 * c5(0) - 32 * c8(0) + 192 * c9(0) ) +
559 c8(1) * ( 12 * c1(0) - 4 * c2(0) + 12 * c3(0) + 108 * c4(0) + 32 * c5(0) + 32 * c6(0) - 192 * c9(0) ) +
560 c9(1) * ( -16 * c1(0) + 16 * c2(0) - 144 * c3(0) + 144 * c4(0) - 192 * c6(0) + 192 * c8(0) )
561 ) +
562 c8(2) * (
563 c1(1) * ( 18 * c2(0) - 2 * c3(0) - 108 * c5(0) + 24 * c6(0) + 12 * c7(0) - 144 * c9(0) ) +
564 c2(1) * ( -18 * c1(0) - 4 * c3(0) - 2 * c4(0) - 12 * c5(0) + 24 * c6(0) - 4 * c7(0) + 16 * c9(0) ) +
565 c3(1) * ( 2 * c1(0) + 4 * c2(0) + 18 * c4(0) + 4 * c5(0) - 24 * c6(0) + 12 * c7(0) - 16 * c9(0) ) +
566 c4(1) * ( 2 * c2(0) - 18 * c3(0) - 12 * c5(0) - 24 * c6(0) + 108 * c7(0) + 144 * c9(0) ) +
567 c5(1) * ( 108 * c1(0) + 12 * c2(0) - 4 * c3(0) + 12 * c4(0) + 32 * c6(0) + 32 * c7(0) - 192 * c9(0) ) +
568 c6(1) * ( -24 * c1(0) - 24 * c2(0) + 24 * c3(0) + 24 * c4(0) - 32 * c5(0) + 32 * c7(0) ) +
569 c7(1) * ( -12 * c1(0) + 4 * c2(0) - 12 * c3(0) - 108 * c4(0) - 32 * c5(0) - 32 * c6(0) + 192 * c9(0) ) +
570 c9(1) * ( 144 * c1(0) - 16 * c2(0) + 16 * c3(0) - 144 * c4(0) + 192 * c5(0) - 192 * c7(0) )
571 ) +
572 c9(2) * (
573 c1(1) * ( 24 * c2(0) - 24 * c4(0) - 144 * c5(0) + 16 * c6(0) - 16 * c7(0) + 144 * c8(0) ) +
574 c2(1) * ( -24 * c1(0) + 24 * c3(0) + 144 * c5(0) - 144 * c6(0) + 16 * c7(0) - 16 * c8(0) ) +
575 c3(1) * ( -24 * c2(0) + 24 * c4(0) - 16 * c5(0) + 144 * c6(0) - 144 * c7(0) + 16 * c8(0) ) +
576 c4(1) * ( 24 * c1(0) - 24 * c3(0) + 16 * c5(0) - 16 * c6(0) + 144 * c7(0) - 144 * c8(0) ) +
577 c5(1) * ( 144 * c1(0) - 144 * c2(0) + 16 * c3(0) - 16 * c4(0) - 192 * c6(0) + 192 * c8(0) ) +
578 c6(1) * ( -16 * c1(0) + 144 * c2(0) - 144 * c3(0) + 16 * c4(0) + 192 * c5(0) - 192 * c7(0) ) +
579 c7(1) * ( 16 * c1(0) - 16 * c2(0) + 144 * c3(0) - 144 * c4(0) + 192 * c6(0) - 192 * c8(0) ) +
580 c8(1) * ( -144 * c1(0) + 16 * c2(0) - 16 * c3(0) + 144 * c4(0) - 192 * c5(0) + 192 * c7(0) ) )
581 ) / 300.0;
582}
583
584std::unique_ptr<IntegrationRule>
585FEI3dHexaTriQuad :: giveIntegrationRule(int order, Element_Geometry_Type egt) const
586{
587 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
589 int points = iRule->getRequiredNumberOfIntegrationPoints(_Cube, order + 15);
590 iRule->SetUpPointsOnCube(points, _Unknown);
591 return std::move(iRule);
592}
593
594std::unique_ptr<IntegrationRule>
595FEI3dHexaTriQuad :: giveBoundaryIntegrationRule(int order, int boundary, Element_Geometry_Type egt) const
596{
597 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
599 int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 6);
600 iRule->SetUpPointsOnSquare(points, _Unknown);
601 return std::move(iRule);
602}
603} // end namespace oofem
oofem::FEI3dHexaTriQuad::evaldNdx
static std::pair< double, FloatMatrixF< 3, 27 > > evaldNdx(const FloatArrayF< 3 > &lcoords, const FEICellGeometry &cellgeo)
Definition fei3dhexatriquad.C:455
oofem::FEI3dHexaTriQuad::computeLocalSurfaceMapping
IntArray computeLocalSurfaceMapping(int iSurf) const override
Definition fei3dhexatriquad.C:237
oofem::FEI3dHexaTriQuad::evalN
static FloatArrayF< 27 > evalN(const FloatArrayF< 3 > &lcoords)
Definition fei3dhexatriquad.C:46
oofem::FEI3dHexaTriQuad::evaldNdxi
static FloatMatrixF< 3, 27 > evaldNdxi(const FloatArrayF< 3 > &lcoords)
Definition fei3dhexatriquad.C:258
oofem::FEICellGeometry::giveVertexCoordinates
virtual const FloatArray giveVertexCoordinates(int i) const =0
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::normalize_giveNorm
double normalize_giveNorm()
Definition floatarray.C:850
oofem::FloatArray::beVectorProductOf
void beVectorProductOf(const FloatArray &v1, const FloatArray &v2)
Definition floatarray.C:476
oofem::FloatArray::add
void add(const FloatArray &src)
Definition floatarray.C:218
oofem::FloatMatrixF::setColumn
void setColumn(const FloatArrayF< N > &src, std::size_t c)
Definition floatmatrixf.h:274
oofem::FloatMatrix::resize
void resize(Index rows, Index cols)
Definition floatmatrix.C:79
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::transpose
FloatMatrixF< M, N > transpose(const FloatMatrixF< N, M > &mat)
Constructs transposed matrix.
Definition floatmatrixf.h:628
oofem::dotT
FloatMatrixF< N, P > dotT(const FloatMatrixF< N, M > &a, const FloatMatrixF< P, M > &b)
Computes .
Definition floatmatrixf.h:659
oofem::det
double det(const FloatMatrixF< 2, 2 > &mat)
Computes the determinant.
Definition floatmatrixf.h:1092
oofem::dot
double dot(const FloatArray &x, const FloatArray &y)
Definition floatarray.C:1011
oofem::inv
FloatMatrixF< N, N > inv(const FloatMatrixF< N, N > &mat, double zeropiv=1e-24)
Computes the inverse.
Definition floatmatrixf.h:1138
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011