OOFEM 3.0
Loading...
Searching...
No Matches
fei3dhexalin.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 "fei3dhexalin.h"
36#include "mathfem.h"
37#include "floatmatrix.h"
38#include "floatarray.h"
39#include "gaussintegrationrule.h"
40#include "floatarrayf.h"
41#include "floatmatrixf.h"
42
43namespace oofem {
44
45FloatArrayF<8>
46FEI3dHexaLin :: evalN(const FloatArrayF<3> &lcoords)
47{
49 double x = lcoords[0];
50 double y = lcoords[1];
51 double z = lcoords[2];
52
53 return {
54 0.125 * ( 1. - x ) * ( 1. - y ) * ( 1. + z ),
55 0.125 * ( 1. - x ) * ( 1. + y ) * ( 1. + z ),
56 0.125 * ( 1. + x ) * ( 1. + y ) * ( 1. + z ),
57 0.125 * ( 1. + x ) * ( 1. - y ) * ( 1. + z ),
58 0.125 * ( 1. - x ) * ( 1. - y ) * ( 1. - z ),
59 0.125 * ( 1. - x ) * ( 1. + y ) * ( 1. - z ),
60 0.125 * ( 1. + x ) * ( 1. + y ) * ( 1. - z ),
61 0.125 * ( 1. + x ) * ( 1. - y ) * ( 1. - z )
62 };
63}
64
65
66void
67FEI3dHexaLin :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
68{
69#if 0
70 answer = evalN(lcoords);
71#else
72 double x = lcoords[0];
73 double y = lcoords[1];
74 double z = lcoords[2];
75
76 answer = Vec8(
77 0.125 * ( 1. - x ) * ( 1. - y ) * ( 1. + z ),
78 0.125 * ( 1. - x ) * ( 1. + y ) * ( 1. + z ),
79 0.125 * ( 1. + x ) * ( 1. + y ) * ( 1. + z ),
80 0.125 * ( 1. + x ) * ( 1. - y ) * ( 1. + z ),
81 0.125 * ( 1. - x ) * ( 1. - y ) * ( 1. - z ),
82 0.125 * ( 1. - x ) * ( 1. + y ) * ( 1. - z ),
83 0.125 * ( 1. + x ) * ( 1. + y ) * ( 1. - z ),
84 0.125 * ( 1. + x ) * ( 1. - y ) * ( 1. - z )
85 );
86#endif
87}
88
89
90FloatMatrixF<3,8>
91FEI3dHexaLin :: evaldNdxi(const FloatArrayF<3> &lcoords)
92{
93 //auto [u, v, w] = lcoords;
94 double u = lcoords.at(1);
95 double v = lcoords.at(2);
96 double w = lcoords.at(3);
97
98 return {
99 -0.125 * ( 1. - v ) * ( 1. + w ),
100 -0.125 * ( 1. - u ) * ( 1. + w ),
101 0.125 * ( 1. - u ) * ( 1. - v ),
102 -0.125 * ( 1. + v ) * ( 1. + w ),
103 0.125 * ( 1. - u ) * ( 1. + w ),
104 0.125 * ( 1. - u ) * ( 1. + v ),
105 0.125 * ( 1. + v ) * ( 1. + w ),
106 0.125 * ( 1. + u ) * ( 1. + w ),
107 0.125 * ( 1. + u ) * ( 1. + v ),
108 0.125 * ( 1. - v ) * ( 1. + w ),
109 -0.125 * ( 1. + u ) * ( 1. + w ),
110 0.125 * ( 1. + u ) * ( 1. - v ),
111 -0.125 * ( 1. - v ) * ( 1. - w ),
112 -0.125 * ( 1. - u ) * ( 1. - w ),
113 -0.125 * ( 1. - u ) * ( 1. - v ),
114 -0.125 * ( 1. + v ) * ( 1. - w ),
115 0.125 * ( 1. - u ) * ( 1. - w ),
116 -0.125 * ( 1. - u ) * ( 1. + v ),
117 0.125 * ( 1. + v ) * ( 1. - w ),
118 0.125 * ( 1. + u ) * ( 1. - w ),
119 -0.125 * ( 1. + u ) * ( 1. + v ),
120 0.125 * ( 1. - v ) * ( 1. - w ),
121 -0.125 * ( 1. + u ) * ( 1. - w ),
122 -0.125 * ( 1. + u ) * ( 1. - v ),
123 };
124}
125
126void
127FEI3dHexaLin :: evaldNdxi(FloatMatrix &dN, const FloatArray &lcoords, const FEICellGeometry &) const
128{
129#if 0
130 dN = evaldNdxi(lcoords);
131#else
132 double u, v, w;
133 u = lcoords.at(1);
134 v = lcoords.at(2);
135 w = lcoords.at(3);
136
137 dN.resize(8, 3);
138
139 dN.at(1, 1) = -0.125 * ( 1. - v ) * ( 1. + w );
140 dN.at(2, 1) = -0.125 * ( 1. + v ) * ( 1. + w );
141 dN.at(3, 1) = 0.125 * ( 1. + v ) * ( 1. + w );
142 dN.at(4, 1) = 0.125 * ( 1. - v ) * ( 1. + w );
143 dN.at(5, 1) = -0.125 * ( 1. - v ) * ( 1. - w );
144 dN.at(6, 1) = -0.125 * ( 1. + v ) * ( 1. - w );
145 dN.at(7, 1) = 0.125 * ( 1. + v ) * ( 1. - w );
146 dN.at(8, 1) = 0.125 * ( 1. - v ) * ( 1. - w );
147
148 dN.at(1, 2) = -0.125 * ( 1. - u ) * ( 1. + w );
149 dN.at(2, 2) = 0.125 * ( 1. - u ) * ( 1. + w );
150 dN.at(3, 2) = 0.125 * ( 1. + u ) * ( 1. + w );
151 dN.at(4, 2) = -0.125 * ( 1. + u ) * ( 1. + w );
152 dN.at(5, 2) = -0.125 * ( 1. - u ) * ( 1. - w );
153 dN.at(6, 2) = 0.125 * ( 1. - u ) * ( 1. - w );
154 dN.at(7, 2) = 0.125 * ( 1. + u ) * ( 1. - w );
155 dN.at(8, 2) = -0.125 * ( 1. + u ) * ( 1. - w );
156
157 dN.at(1, 3) = 0.125 * ( 1. - u ) * ( 1. - v );
158 dN.at(2, 3) = 0.125 * ( 1. - u ) * ( 1. + v );
159 dN.at(3, 3) = 0.125 * ( 1. + u ) * ( 1. + v );
160 dN.at(4, 3) = 0.125 * ( 1. + u ) * ( 1. - v );
161 dN.at(5, 3) = -0.125 * ( 1. - u ) * ( 1. - v );
162 dN.at(6, 3) = -0.125 * ( 1. - u ) * ( 1. + v );
163 dN.at(7, 3) = -0.125 * ( 1. + u ) * ( 1. + v );
164 dN.at(8, 3) = -0.125 * ( 1. + u ) * ( 1. - v );
165#endif
166}
167
168
169std::pair<double, FloatMatrixF<3,8>>
170FEI3dHexaLin :: evaldNdx(const FloatArrayF<3> &lcoords, const FEICellGeometry &cellgeo)
171{
172 auto dNduvw = evaldNdxi(lcoords);
173 FloatMatrixF<3,8> coords;
174 for ( int i = 0; i < 8; i++ ) {
176 coords.setColumn(cellgeo.giveVertexCoordinates(i+1), i);
177 }
178 auto jacT = dotT(dNduvw, coords);
179 return {det(jacT), dot(inv(jacT), dNduvw)};
180}
181
182
183double
184FEI3dHexaLin :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
185{
186#if 0
187 auto tmp = evaldNdx(lcoords, cellgeo);
188 //auto [det, dndx] = evaldNdx(lcoords, cellgeo);
189 answer = tmp.second;
190 return tmp.first;
191#else
192 FloatMatrix jacobianMatrix, inv, dNduvw, coords;
193
194 this->evaldNdxi(dNduvw, lcoords, cellgeo);
195 coords.resize(3, 8);
196 for ( int i = 1; i <= 8; i++ ) {
197 coords.setColumn(cellgeo.giveVertexCoordinates(i), i);
198 }
199 jacobianMatrix.beProductOf(coords, dNduvw);
200 inv.beInverseOf(jacobianMatrix);
201
202 answer.beProductOf(dNduvw, inv);
203 return jacobianMatrix.giveDeterminant();
204#endif
205}
206
207void
208FEI3dHexaLin :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
209{
210 FloatArray n;
211 this->evalN(n, lcoords, cellgeo);
212
213 answer.clear();
214 for ( int i = 1; i <= 8; i++ ) {
215 answer.add( n.at(i), cellgeo.giveVertexCoordinates(i) );
216 }
217}
218
219#define POINT_TOL 1.e-3
220
221int
222FEI3dHexaLin :: global2local(FloatArray &answer, const FloatArray &coords, const FEICellGeometry &cellgeo) const
223{
224 double x1, x2, x3, x4, x5, x6, x7, x8, a1, a2, a3, a4, a5, a6, a7, a8;
225 double y1, y2, y3, y4, y5, y6, y7, y8, b1, b2, b3, b4, b5, b6, b7, b8;
226 double z1, z2, z3, z4, z5, z6, z7, z8, c1, c2, c3, c4, c5, c6, c7, c8;
227 double xp, yp, zp, u, v, w;
228 FloatMatrix p(3, 3);
229 FloatArray r(3), delta;
230 int nite = 0;
231
232 x1 = cellgeo.giveVertexCoordinates(1).at(1);
233 x2 = cellgeo.giveVertexCoordinates(2).at(1);
234 x3 = cellgeo.giveVertexCoordinates(3).at(1);
235 x4 = cellgeo.giveVertexCoordinates(4).at(1);
236 x5 = cellgeo.giveVertexCoordinates(5).at(1);
237 x6 = cellgeo.giveVertexCoordinates(6).at(1);
238 x7 = cellgeo.giveVertexCoordinates(7).at(1);
239 x8 = cellgeo.giveVertexCoordinates(8).at(1);
240
241 y1 = cellgeo.giveVertexCoordinates(1).at(2);
242 y2 = cellgeo.giveVertexCoordinates(2).at(2);
243 y3 = cellgeo.giveVertexCoordinates(3).at(2);
244 y4 = cellgeo.giveVertexCoordinates(4).at(2);
245 y5 = cellgeo.giveVertexCoordinates(5).at(2);
246 y6 = cellgeo.giveVertexCoordinates(6).at(2);
247 y7 = cellgeo.giveVertexCoordinates(7).at(2);
248 y8 = cellgeo.giveVertexCoordinates(8).at(2);
249
250 z1 = cellgeo.giveVertexCoordinates(1).at(3);
251 z2 = cellgeo.giveVertexCoordinates(2).at(3);
252 z3 = cellgeo.giveVertexCoordinates(3).at(3);
253 z4 = cellgeo.giveVertexCoordinates(4).at(3);
254 z5 = cellgeo.giveVertexCoordinates(5).at(3);
255 z6 = cellgeo.giveVertexCoordinates(6).at(3);
256 z7 = cellgeo.giveVertexCoordinates(7).at(3);
257 z8 = cellgeo.giveVertexCoordinates(8).at(3);
258
259 xp = coords.at(1);
260 yp = coords.at(2);
261 zp = coords.at(3);
262
263 a1 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8;
264 a2 = -x1 - x2 + x3 + x4 - x5 - x6 + x7 + x8;
265 a3 = -x1 + x2 + x3 - x4 - x5 + x6 + x7 - x8;
266 a4 = x1 + x2 + x3 + x4 - x5 - x6 - x7 - x8;
267 a5 = x1 - x2 + x3 - x4 + x5 - x6 + x7 - x8;
268 a6 = -x1 - x2 + x3 + x4 + x5 + x6 - x7 - x8;
269 a7 = -x1 + x2 + x3 - x4 + x5 - x6 - x7 + x8;
270 a8 = x1 - x2 + x3 - x4 - x5 + x6 - x7 + x8;
271
272 b1 = y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8;
273 b2 = -y1 - y2 + y3 + y4 - y5 - y6 + y7 + y8;
274 b3 = -y1 + y2 + y3 - y4 - y5 + y6 + y7 - y8;
275 b4 = y1 + y2 + y3 + y4 - y5 - y6 - y7 - y8;
276 b5 = y1 - y2 + y3 - y4 + y5 - y6 + y7 - y8;
277 b6 = -y1 - y2 + y3 + y4 + y5 + y6 - y7 - y8;
278 b7 = -y1 + y2 + y3 - y4 + y5 - y6 - y7 + y8;
279 b8 = y1 - y2 + y3 - y4 - y5 + y6 - y7 + y8;
280
281 c1 = z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8;
282 c2 = -z1 - z2 + z3 + z4 - z5 - z6 + z7 + z8;
283 c3 = -z1 + z2 + z3 - z4 - z5 + z6 + z7 - z8;
284 c4 = z1 + z2 + z3 + z4 - z5 - z6 - z7 - z8;
285 c5 = z1 - z2 + z3 - z4 + z5 - z6 + z7 - z8;
286 c6 = -z1 - z2 + z3 + z4 + z5 + z6 - z7 - z8;
287 c7 = -z1 + z2 + z3 - z4 + z5 - z6 - z7 + z8;
288 c8 = z1 - z2 + z3 - z4 - z5 + z6 - z7 + z8;
289
290 // setup initial guess
291 answer.resize(3);
292 answer.zero();
293
294 // apply Newton-Raphson to solve the problem
295 for ( ;; ) {
296 if ( ( ++nite ) > 10 ) {
297 answer.zero();
298 return false;
299 }
300
301 u = answer.at(1);
302 v = answer.at(2);
303 w = answer.at(3);
304
305 // compute the residual
306 r.at(1) = a1 + u * a2 + v * a3 + w * a4 + u * v * a5 + u * w * a6 + v * w * a7 + u * v * w * a8 - 8.0 * xp;
307 r.at(2) = b1 + u * b2 + v * b3 + w * b4 + u * v * b5 + u * w * b6 + v * w * b7 + u * v * w * b8 - 8.0 * yp;
308 r.at(3) = c1 + u * c2 + v * c3 + w * c4 + u * v * c5 + u * w * c6 + v * w * c7 + u * v * w * c8 - 8.0 * zp;
309
310 // check for convergence
311 if ( r.computeSquaredNorm() < 1.e-20 ) {
312 break; // sqrt(1.e-20) = 1.e-10
313 }
314
315 p.at(1, 1) = a2 + v * a5 + w * a6 + v * w * a8;
316 p.at(1, 2) = a3 + u * a5 + w * a7 + u * w * a8;
317 p.at(1, 3) = a4 + u * a6 + v * a7 + u * v * a8;
318
319 p.at(2, 1) = b2 + v * b5 + w * b6 + v * w * b8;
320 p.at(2, 2) = b3 + u * b5 + w * b7 + u * w * b8;
321 p.at(2, 3) = b4 + u * b6 + v * b7 + u * v * b8;
322
323 p.at(3, 1) = c2 + v * c5 + w * c6 + v * w * c8;
324 p.at(3, 2) = c3 + u * c5 + w * c7 + u * w * c8;
325 p.at(3, 3) = c4 + u * c6 + v * c7 + u * v * c8;
326
327 // solve for corrections
328 p.solveForRhs(r, delta);
329
330 // update guess
331 answer.subtract(delta);
332 }
333
334 // test if inside
335 bool inside = true;
336 for ( int i = 1; i <= 3; i++ ) {
337 if ( answer.at(i) < ( -1. - POINT_TOL ) ) {
338 answer.at(i) = -1.;
339 inside = false;
340 } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
341 answer.at(i) = 1.;
342 inside = false;
343 }
344 }
345
346 return inside;
347}
348
349void
350FEI3dHexaLin :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
351{
352 double ksi = lcoords.at(1);
353 answer.resize(2);
354
355 answer.at(1) = ( 1. - ksi ) * 0.5;
356 answer.at(2) = ( 1. + ksi ) * 0.5;
357}
358
359void
360FEI3dHexaLin :: edgeEvaldNdx(FloatMatrix &answer, int iedge,
361 const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
362{
363 const auto &edgeNodes = this->computeLocalEdgeMapping(iedge);
364 double l = this->edgeComputeLength(edgeNodes, cellgeo);
365
366 answer.resize(2, 1);
367 answer.at(1, 1) = -1.0 / l;
368 answer.at(2, 1) = 1.0 / l;
369}
370
371
372void
373FEI3dHexaLin :: edgeEvaldNdxi(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
374{
375 answer.resize(2);
376 answer(0) = -0.5;
377 answer(1) = 0.5;
378}
379
380
381void
382FEI3dHexaLin :: edgeLocal2global(FloatArray &answer, int iedge,
383 const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
384{
385 FloatArray n;
386 const auto &edgeNodes = this->computeLocalEdgeMapping(iedge);
387 this->edgeEvalN(n, iedge, lcoords, cellgeo);
388
389 answer.resize(3);
390 answer.at(1) = n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) ).at(1) +
391 n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) ).at(1);
392 answer.at(2) = n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) ).at(2) +
393 n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) ).at(2);
394 answer.at(3) = n.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) ).at(3) +
395 n.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) ).at(3);
396}
397
398
399double
400FEI3dHexaLin :: edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
401{
402 const auto &edgeNodes = this->computeLocalEdgeMapping(iedge);
403 return 0.5 * this->edgeComputeLength(edgeNodes, cellgeo);
404}
405
406
407IntArray
408FEI3dHexaLin :: computeLocalEdgeMapping(int iedge) const
409{
410 if ( iedge == 1 ) { // edge between nodes 1 2
411 return {1, 2};
412 } else if ( iedge == 2 ) { // edge between nodes 2 3
413 return {2, 3};
414 } else if ( iedge == 3 ) { // edge between nodes 3 4
415 return {3, 4};
416 } else if ( iedge == 4 ) { // edge between nodes 4 1
417 return {4, 1};
418 } else if ( iedge == 5 ) { // edge between nodes 1 5
419 return {1, 5};
420 } else if ( iedge == 6 ) { // edge between nodes 2 6
421 return {2, 6};
422 } else if ( iedge == 7 ) { // edge between nodes 3 7
423 return {3, 7};
424 } else if ( iedge == 8 ) { // edge between nodes 4 8
425 return {4, 8};
426 } else if ( iedge == 9 ) { // edge between nodes 5 6
427 return {5, 6};
428 } else if ( iedge == 10 ) { // edge between nodes 6 7
429 return {6, 7};
430 } else if ( iedge == 11 ) { // edge between nodes 7 8
431 return {7, 8};
432 } else if ( iedge == 12 ) { // edge between nodes 8 5
433 return {8, 5};
434 } else {
435 throw std::range_error("invalid edge number");
436 //return {};
437 }
438}
439
440double
441FEI3dHexaLin :: edgeComputeLength(const IntArray &edgeNodes, const FEICellGeometry &cellgeo) const
442{
443 return distance(cellgeo.giveVertexCoordinates( edgeNodes.at(2) ), cellgeo.giveVertexCoordinates( edgeNodes.at(1) ));
444}
445
446void
447FEI3dHexaLin :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
448{
449 double ksi = lcoords.at(1);
450 double eta = lcoords.at(2);
451
452 answer.resize(4);
453 answer.at(1) = ( 1. + ksi ) * ( 1. + eta ) * 0.25;
454 answer.at(2) = ( 1. - ksi ) * ( 1. + eta ) * 0.25;
455 answer.at(3) = ( 1. - ksi ) * ( 1. - eta ) * 0.25;
456 answer.at(4) = ( 1. + ksi ) * ( 1. - eta ) * 0.25;
457}
458
459void
460FEI3dHexaLin :: surfaceEvaldNdx(FloatMatrix &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
461{
462 // Note, this must be in correct order, not just the correct nodes, therefore we must use snodes;
463 const auto &snodes = this->computeLocalSurfaceMapping(isurf);
464
465 FloatArray lcoords_hex;
466
468#if 1
469 if ( isurf == 1 ) { // surface 1 - nodes 1 4 3 2
470 lcoords_hex = Vec3(-lcoords.at(1), -lcoords.at(2), 1);
471 } else if ( isurf == 2 ) { // surface 2 - nodes 5 6 7 8
472 lcoords_hex = Vec3(-lcoords.at(2), -lcoords.at(1), -1);
473 } else if ( isurf == 3 ) { // surface 3 - nodes 1 2 6 5
474 lcoords_hex = Vec3(-1, -lcoords.at(1), lcoords.at(2));
475 } else if ( isurf == 4 ) { // surface 4 - nodes 2 3 7 6
476 lcoords_hex = Vec3(-lcoords.at(1), 1, lcoords.at(2));
477 } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7
478 lcoords_hex = Vec3(1, lcoords.at(1), lcoords.at(2));
479 } else if ( isurf == 6 ) { // surface 6 - nodes 4 1 5 8
480 lcoords_hex = Vec3(lcoords.at(1), -1, lcoords.at(2));
481 } else {
482 OOFEM_ERROR("wrong surface number (%d)", isurf);
483 }
484#else
486 if ( isurf == 1 ) { // surface 1 - nodes 3 4 8 7
487 lcoords_hex = Vec3(-1, lcoords.at(1), lcoords.at(2));
488 } else if ( isurf == 2 ) { // surface 2 - nodes 2 1 5 6
489 lcoords_hex = Vec3(1, lcoords.at(1), lcoords.at(2));
490 } else if ( isurf == 3 ) { // surface 3 - nodes 3 7 6 2
491 lcoords_hex = Vec3(lcoords.at(1), -1, lcoords.at(2));
492 } else if ( isurf == 4 ) { // surface 4 - nodes 4 8 5 1
493 lcoords_hex = Vec3(lcoords.at(1), 1, lcoords.at(2));
494 } else if ( isurf == 5 ) { // surface 5 - nodes 3 2 1 4
495 lcoords_hex = Vec3(lcoords.at(1), lcoords.at(2), -1);
496 } else if ( isurf == 6 ) { // surface 6 - nodes 7 6 5 8
497 lcoords_hex = Vec3(lcoords.at(1), lcoords.at(2), 1);
498 } else {
499 OOFEM_ERROR("wrong surface number (%d)", isurf);
500 }
501#endif
502
503 FloatMatrix fullB;
504 this->evaldNdx(fullB, lcoords_hex, cellgeo);
505 answer.resize(snodes.giveSize(), 3);
506 for ( int i = 1; i <= snodes.giveSize(); ++i ) {
507 for ( int j = 1; j <= 3; ++j ) {
508 answer.at(i, j) = fullB.at(snodes.at(i), j);
509 }
510 }
511}
512
513void
514FEI3dHexaLin :: surfaceLocal2global(FloatArray &answer, int iedge,
515 const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
516{
517 FloatArray n;
518
519 const auto &nodes = this->computeLocalSurfaceMapping(iedge);
520 this->surfaceEvalN(n, iedge, lcoords, cellgeo);
521
522 answer.resize(3);
523 answer.at(1) = n.at(1) * cellgeo.giveVertexCoordinates( nodes.at(1) ).at(1) + n.at(2) * cellgeo.giveVertexCoordinates( nodes.at(2) ).at(1) +
524 n.at(3) * cellgeo.giveVertexCoordinates( nodes.at(3) ).at(1) + n.at(4) * cellgeo.giveVertexCoordinates( nodes.at(4) ).at(1);
525 answer.at(2) = n.at(1) * cellgeo.giveVertexCoordinates( nodes.at(1) ).at(2) + n.at(2) * cellgeo.giveVertexCoordinates( nodes.at(2) ).at(2) +
526 n.at(3) * cellgeo.giveVertexCoordinates( nodes.at(3) ).at(2) + n.at(4) * cellgeo.giveVertexCoordinates( nodes.at(4) ).at(2);
527 answer.at(3) = n.at(1) * cellgeo.giveVertexCoordinates( nodes.at(1) ).at(3) + n.at(2) * cellgeo.giveVertexCoordinates( nodes.at(2) ).at(3) +
528 n.at(3) * cellgeo.giveVertexCoordinates( nodes.at(3) ).at(3) + n.at(4) * cellgeo.giveVertexCoordinates( nodes.at(4) ).at(3);
529}
530
531double
532FEI3dHexaLin :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
533{
534 FloatArray a, b, dNdksi(4), dNdeta(4);
535 const auto &snodes = this->computeLocalSurfaceMapping(isurf);
536
537 double ksi = lcoords.at(1);
538 double eta = lcoords.at(2);
539
540 // No need to divide by 1/4, we'll normalize anyway;
541 dNdksi.at(1) = ( 1. + eta );
542 dNdksi.at(2) = -( 1. + eta );
543 dNdksi.at(3) = -( 1. - eta );
544 dNdksi.at(4) = ( 1. - eta );
545
546 dNdeta.at(1) = ( 1. + ksi );
547 dNdeta.at(2) = ( 1. - ksi );
548 dNdeta.at(3) = -( 1. - ksi );
549 dNdeta.at(4) = -( 1. + ksi );
550
551 for ( int i = 1; i <= 4; ++i ) {
552 a.add( dNdksi.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
553 b.add( dNdeta.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
554 }
555
556 answer.beVectorProductOf(a, b);
557 return answer.normalize_giveNorm() * 0.0625;
558}
559
560double
561FEI3dHexaLin :: surfaceGiveTransformationJacobian(int isurf, const FloatArray &lcoords,
562 const FEICellGeometry &cellgeo) const
563{
564 FloatArray normal;
565 return this->surfaceEvalNormal(normal, isurf, lcoords, cellgeo);
566}
567
568void
569FEI3dHexaLin :: surfaceEvaldNdxi(FloatMatrix &answer, const FloatArray &lcoords) const
570{
571 // Returns matrix with derivatives wrt local coordinates
572 answer.resize(4, 2);
573 double ksi = lcoords.at(1);
574 double eta = lcoords.at(2);
575
576 for ( int i = 1; i <= 3; ++i ) {
577 answer.at(1,1) = 0.25 * ( 1. + eta );
578 answer.at(2,1) = -0.25 * ( 1. + eta );
579 answer.at(3,1) = -0.25 * ( 1. - eta );
580 answer.at(4,1) = 0.25 * ( 1. - eta );
581
582 answer.at(1,2) = 0.25 * ( 1. + ksi );
583 answer.at(2,2) = 0.25 * ( 1. - ksi );
584 answer.at(3,2) = -0.25 * ( 1. - ksi );
585 answer.at(4,2) = -0.25 * ( 1. + ksi );
586 }
587}
588
589
590
591
592
593IntArray
594FEI3dHexaLin :: computeLocalSurfaceMapping(int isurf) const
595{
596 if ( isurf == 1 ) { // surface 1 - nodes 1 4 3 2
597 return {1, 4, 3, 2};
598 } else if ( isurf == 2 ) { // surface 2 - nodes 5 6 7 8
599 return {5, 6, 7, 8};
600 } else if ( isurf == 3 ) { // surface 3 - nodes 1 2 6 5
601 return {1, 2, 6, 5};
602 } else if ( isurf == 4 ) { // surface 4 - nodes 2 3 7 6
603 return {2, 3, 7, 6};
604 } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7
605 return {3, 4, 8, 7};
606 } else if ( isurf == 6 ) { // surface 6 - nodes 4 1 5 8
607 return {4, 1, 5, 8};
608 } else {
609 throw std::runtime_error("invalid surface number");
610 }
611}
612
613void
614FEI3dHexaLin :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
615// Returns the jacobian matrix J (x,y,z)/(ksi,eta,dzeta) of the receiver.
616{
617 FloatMatrix dNduvw, coords;
618 this->evaldNdxi(dNduvw, lcoords, cellgeo);
619 coords.resize(3, 8);
620 for ( int i = 1; i <= 8; i++ ) {
621 coords.setColumn(cellgeo.giveVertexCoordinates(i), i);
622 }
623 jacobianMatrix.beProductOf(coords, dNduvw);
624}
625
626
627double
628FEI3dHexaLin :: evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo) const
629{
630 const auto &fNodes = this->computeLocalSurfaceMapping(iEdge);
631
632 const auto &c1 = cellgeo.giveVertexCoordinates( fNodes.at(1) );
633 const auto &c2 = cellgeo.giveVertexCoordinates( fNodes.at(2) );
634 const auto &c3 = cellgeo.giveVertexCoordinates( fNodes.at(3) );
635 const auto &c4 = cellgeo.giveVertexCoordinates( fNodes.at(4) );
636
637 return (
638 c4(2) * ( c1(1) * ( -c2(0) - c3(0) ) + c2(1) * ( c1(0) - c3(0) ) + c3(1) * ( c1(0) + c2(0) ) ) +
639 c3(2) * ( c1(1) * ( -c2(0) + c4(0) ) + c2(1) * ( c1(0) + c4(0) ) + c4(1) * ( -c1(0) - c2(0) ) ) +
640 c2(2) * ( c1(1) * ( c3(0) + c4(0) ) + c3(1) * ( -c1(0) - c4(0) ) + c4(1) * ( -c1(0) + c3(0) ) ) +
641 c1(2) * ( c2(1) * ( -c3(0) - c4(0) ) + c3(1) * ( c2(0) - c4(0) ) + c4(1) * ( c2(0) + c3(0) ) ) ) * 0.25;
642}
643
644std::unique_ptr<IntegrationRule>
645FEI3dHexaLin :: giveIntegrationRule(int order, Element_Geometry_Type egt) const
646{
647 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
648 int points = iRule->getRequiredNumberOfIntegrationPoints(_Cube, order + 6);
649 iRule->SetUpPointsOnCube(points, _Unknown);
650 return std::move(iRule);
651}
652
653std::unique_ptr<IntegrationRule>
654FEI3dHexaLin :: giveBoundaryIntegrationRule(int order, int boundary, Element_Geometry_Type egt) const
655{
656 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
657 int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 2);
658 iRule->SetUpPointsOnSquare(points, _Unknown);
659 return std::move(iRule);
660}
661} // end namespace oofem
oofem::FEI3dHexaLin::edgeComputeLength
double edgeComputeLength(const IntArray &edgeNodes, const FEICellGeometry &cellgeo) const
Definition fei3dhexalin.C:441
oofem::FEI3dHexaLin::evaldNdxi
static FloatMatrixF< 3, 8 > evaldNdxi(const FloatArrayF< 3 > &lcoords)
Definition fei3dhexalin.C:91
oofem::FEI3dHexaLin::computeLocalSurfaceMapping
IntArray computeLocalSurfaceMapping(int iedge) const override
Definition fei3dhexalin.C:594
oofem::FEI3dHexaLin::evaldNdx
static std::pair< double, FloatMatrixF< 3, 8 > > evaldNdx(const FloatArrayF< 3 > &lcoords, const FEICellGeometry &cellgeo)
Definition fei3dhexalin.C:170
oofem::FEI3dHexaLin::edgeEvalN
void edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dhexalin.C:350
oofem::FEI3dHexaLin::surfaceEvalN
void surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dhexalin.C:447
oofem::FEI3dHexaLin::evalN
static FloatArrayF< 8 > evalN(const FloatArrayF< 3 > &lcoords)
Definition fei3dhexalin.C:46
oofem::FEI3dHexaLin::surfaceEvalNormal
double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dhexalin.C:532
oofem::FEI3dHexaLin::computeLocalEdgeMapping
IntArray computeLocalEdgeMapping(int iedge) const override
Definition fei3dhexalin.C:408
oofem::FEICellGeometry::giveVertexCoordinates
virtual const FloatArray giveVertexCoordinates(int i) const =0
oofem::FloatArrayF::at
double & at(std::size_t i)
Definition floatarrayf.h:157
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition floatarray.C:683
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::FloatArray::subtract
void subtract(const FloatArray &src)
Definition floatarray.C:320
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::beProductOf
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Definition floatmatrix.C:365
oofem::FloatMatrix::setColumn
void setColumn(const FloatArray &src, int c)
Definition floatmatrix.C:641
oofem::FloatMatrix::giveDeterminant
double giveDeterminant() const
Definition floatmatrix.C:1085
oofem::FloatMatrix::at
double at(std::size_t i, std::size_t j) const
Definition floatmatrix.h:221
oofem::FloatMatrix::solveForRhs
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Definition floatmatrix.C:1303
oofem::IntArray::at
int & at(std::size_t i)
Definition intarray.h:104
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
POINT_TOL
#define POINT_TOL
Definition fei2dquadlin.C:128
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::Vec8
static FloatArray Vec8(const double &a, const double &b, const double &c, const double &d, const double &e, const double &f, const double &g, const double &h)
Definition floatarray.h:612
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
oofem::distance
double distance(const FloatArray &x, const FloatArray &y)
Definition floatarray.C:1016
oofem::Vec3
static FloatArray Vec3(const double &a, const double &b, const double &c)
Definition floatarray.h:607
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