OOFEM 3.0
Loading...
Searching...
No Matches
fei3dhexaquad.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 "fei3dhexaquad.h"
36#include "intarray.h"
37#include "floatarray.h"
38#include "floatmatrix.h"
39#include "floatarrayf.h"
40#include "floatmatrixf.h"
41#include "gaussintegrationrule.h"
42
43namespace oofem {
44
45FloatArrayF<20>
46FEI3dHexaQuad :: 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 return {
53 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 + w ) * ( -u - v + w - 2.0 ),
54 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 + w ) * ( -u + v + w - 2.0 ),
55 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 + w ) * ( u + v + w - 2.0 ),
56 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 + w ) * ( u - v + w - 2.0 ),
57 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 - w ) * ( -u - v - w - 2.0 ),
58 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 - w ) * ( -u + v - w - 2.0 ),
59 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 - w ) * ( u + v - w - 2.0 ),
60 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 - w ) * ( u - v - w - 2.0 ),
61 0.25 * ( 1.0 - v * v ) * ( 1.0 - u ) * ( 1.0 + w ),
62 0.25 * ( 1.0 - u * u ) * ( 1.0 + v ) * ( 1.0 + w ),
63 0.25 * ( 1.0 - v * v ) * ( 1.0 + u ) * ( 1.0 + w ),
64 0.25 * ( 1.0 - u * u ) * ( 1.0 - v ) * ( 1.0 + w ),
65 0.25 * ( 1.0 - v * v ) * ( 1.0 - u ) * ( 1.0 - w ),
66 0.25 * ( 1.0 - u * u ) * ( 1.0 + v ) * ( 1.0 - w ),
67 0.25 * ( 1.0 - v * v ) * ( 1.0 + u ) * ( 1.0 - w ),
68 0.25 * ( 1.0 - u * u ) * ( 1.0 - v ) * ( 1.0 - w ),
69 0.25 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 - w * w ),
70 0.25 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 - w * w ),
71 0.25 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 - w * w ),
72 0.25 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 - w * w )
73 };
74}
75
76
77void
78FEI3dHexaQuad :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
79{
80#if 0
81 answer = evalN(lcoords);
82#else
83 double u, v, w;
84 answer.resize(20);
85
86 u = lcoords.at(1);
87 v = lcoords.at(2);
88 w = lcoords.at(3);
89
90 answer.at(1) = 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 + w ) * ( -u - v + w - 2.0 );
91 answer.at(2) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 + w ) * ( -u + v + w - 2.0 );
92 answer.at(3) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 + w ) * ( u + v + w - 2.0 );
93 answer.at(4) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 + w ) * ( u - v + w - 2.0 );
94 answer.at(5) = 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 - w ) * ( -u - v - w - 2.0 );
95 answer.at(6) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 - w ) * ( -u + v - w - 2.0 );
96 answer.at(7) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 - w ) * ( u + v - w - 2.0 );
97 answer.at(8) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 - w ) * ( u - v - w - 2.0 );
98
99 answer.at(9) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - u ) * ( 1.0 + w );
100 answer.at(10) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + v ) * ( 1.0 + w );
101 answer.at(11) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + u ) * ( 1.0 + w );
102 answer.at(12) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - v ) * ( 1.0 + w );
103
104 answer.at(13) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - u ) * ( 1.0 - w );
105 answer.at(14) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + v ) * ( 1.0 - w );
106 answer.at(15) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + u ) * ( 1.0 - w );
107 answer.at(16) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - v ) * ( 1.0 - w );
108
109 answer.at(17) = 0.25 * ( 1.0 - u ) * ( 1.0 - v ) * ( 1.0 - w * w );
110 answer.at(18) = 0.25 * ( 1.0 - u ) * ( 1.0 + v ) * ( 1.0 - w * w );
111 answer.at(19) = 0.25 * ( 1.0 + u ) * ( 1.0 + v ) * ( 1.0 - w * w );
112 answer.at(20) = 0.25 * ( 1.0 + u ) * ( 1.0 - v ) * ( 1.0 - w * w );
113#endif
114}
115
116
117FloatMatrixF<3,20>
118FEI3dHexaQuad :: evaldNdxi(const FloatArrayF<3> &lcoords)
119{
120 //auto [u, v, w] = lcoords;
121 double u = lcoords[0];
122 double v = lcoords[1];
123 double w = lcoords[2];
124
125 return {
126 0.125 * ( 1.0 - v ) * ( 1.0 + w ) * ( 2.0 * u + v - w + 1.0 ),
127 0.125 * ( 1.0 - u ) * ( 1.0 + w ) * ( 2.0 * v + u - w + 1.0 ),
128 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 2.0 * w - u - v - 1.0 ),
129 0.125 * ( 1.0 + v ) * ( 1.0 + w ) * ( 2.0 * u - v - w + 1.0 ),
130 0.125 * ( 1.0 - u ) * ( 1.0 + w ) * ( 2.0 * v - u + w - 1.0 ),
131 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 2.0 * w - u + v - 1.0 ),
132 0.125 * ( 1.0 + v ) * ( 1.0 + w ) * ( 2.0 * u + v + w - 1.0 ),
133 0.125 * ( 1.0 + u ) * ( 1.0 + w ) * ( 2.0 * v + u + w - 1.0 ),
134 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 2.0 * w + u + v - 1.0 ),
135 0.125 * ( 1.0 - v ) * ( 1.0 + w ) * ( 2.0 * u - v + w - 1.0 ),
136 0.125 * ( 1.0 + u ) * ( 1.0 + w ) * ( 2.0 * v - u - w + 1.0 ),
137 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 2.0 * w + u - v - 1.0 ),
138 0.125 * ( 1.0 - v ) * ( 1.0 - w ) * ( 2.0 * u + v + w + 1.0 ),
139 0.125 * ( 1.0 - u ) * ( 1.0 - w ) * ( 2.0 * v + u + w + 1.0 ),
140 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 2.0 * w + u + v + 1.0 ),
141 0.125 * ( 1.0 + v ) * ( 1.0 - w ) * ( 2.0 * u - v + w + 1.0 ),
142 0.125 * ( 1.0 - u ) * ( 1.0 - w ) * ( 2.0 * v - u - w - 1.0 ),
143 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 2.0 * w + u - v + 1.0 ),
144 0.125 * ( 1.0 + v ) * ( 1.0 - w ) * ( 2.0 * u + v - w - 1.0 ),
145 0.125 * ( 1.0 + u ) * ( 1.0 - w ) * ( 2.0 * v + u - w - 1.0 ),
146 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 2.0 * w - u - v + 1.0 ),
147 0.125 * ( 1.0 - v ) * ( 1.0 - w ) * ( 2.0 * u - v - w - 1.0 ),
148 0.125 * ( 1.0 + u ) * ( 1.0 - w ) * ( 2.0 * v - u + w + 1.0 ),
149 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 2.0 * w - u + v + 1.0 ),
150 -0.25 * ( 1.0 - v * v ) * ( 1.0 + w ),
151 -0.50 * v * ( 1.0 - u ) * ( 1.0 + w ),
152 0.25 * ( 1.0 - v * v ) * ( 1.0 - u ),
153 -0.50 * u * ( 1.0 + v ) * ( 1.0 + w ),
154 0.25 * ( 1.0 - u * u ) * ( 1.0 + w ),
155 0.25 * ( 1.0 - u * u ) * ( 1.0 + v ),
156 -0.50 * v * ( 1.0 + u ) * ( 1.0 + w ),
157 0.25 * ( 1.0 - v * v ) * ( 1.0 + w ),
158 0.25 * ( 1.0 - v * v ) * ( 1.0 + u ),
159 -0.50 * u * ( 1.0 - v ) * ( 1.0 + w ),
160 -0.25 * ( 1.0 - u * u ) * ( 1.0 + w ),
161 0.25 * ( 1.0 - u * u ) * ( 1.0 - v ),
162 -0.25 * ( 1.0 - v * v ) * ( 1.0 - w ),
163 -0.50 * v * ( 1.0 - u ) * ( 1.0 - w ),
164 -0.25 * ( 1.0 - v * v ) * ( 1.0 - u ),
165 -0.50 * u * ( 1.0 + v ) * ( 1.0 - w ),
166 0.25 * ( 1.0 - u * u ) * ( 1.0 - w ),
167 -0.25 * ( 1.0 - u * u ) * ( 1.0 + v ),
168 0.25 * ( 1.0 - v * v ) * ( 1.0 - w ),
169 -0.50 * v * ( 1.0 + u ) * ( 1.0 - w ),
170 -0.25 * ( 1.0 - v * v ) * ( 1.0 + u ),
171 -0.50 * u * ( 1.0 - v ) * ( 1.0 - w ),
172 -0.25 * ( 1.0 - u * u ) * ( 1.0 - w ),
173 -0.25 * ( 1.0 - u * u ) * ( 1.0 - v ),
174 -0.25 * ( 1.0 - v ) * ( 1.0 - w * w ),
175 -0.25 * ( 1.0 - u ) * ( 1.0 - w * w ),
176 -0.50 * ( 1.0 - u ) * ( 1.0 - v ) * w,
177 -0.25 * ( 1.0 + v ) * ( 1.0 - w * w ),
178 0.25 * ( 1.0 - u ) * ( 1.0 - w * w ),
179 -0.50 * ( 1.0 - u ) * ( 1.0 + v ) * w,
180 0.25 * ( 1.0 + v ) * ( 1.0 - w * w ),
181 0.25 * ( 1.0 + u ) * ( 1.0 - w * w ),
182 -0.50 * ( 1.0 + u ) * ( 1.0 + v ) * w,
183 0.25 * ( 1.0 - v ) * ( 1.0 - w * w ),
184 -0.25 * ( 1.0 + u ) * ( 1.0 - w * w ),
185 -0.50 * ( 1.0 + u ) * ( 1.0 - v ) * w,
186 };
187}
188
189
190void
191FEI3dHexaQuad :: evaldNdxi(FloatMatrix &dN, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
192{
193#if 0
194 dN = evaldNdxi(lcoords);
195#else
196 double u, v, w;
197 u = lcoords.at(1);
198 v = lcoords.at(2);
199 w = lcoords.at(3);
200
201 dN.resize(20, 3);
202
203 dN.at(1, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 + w ) * ( 2.0 * u + v - w + 1.0 );
204 dN.at(2, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 + w ) * ( 2.0 * u - v - w + 1.0 );
205 dN.at(3, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 + w ) * ( 2.0 * u + v + w - 1.0 );
206 dN.at(4, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 + w ) * ( 2.0 * u - v + w - 1.0 );
207 dN.at(5, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 - w ) * ( 2.0 * u + v + w + 1.0 );
208 dN.at(6, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 - w ) * ( 2.0 * u - v + w + 1.0 );
209 dN.at(7, 1) = 0.125 * ( 1.0 + v ) * ( 1.0 - w ) * ( 2.0 * u + v - w - 1.0 );
210 dN.at(8, 1) = 0.125 * ( 1.0 - v ) * ( 1.0 - w ) * ( 2.0 * u - v - w - 1.0 );
211 dN.at(9, 1) = -0.25 * ( 1.0 - v * v ) * ( 1.0 + w );
212 dN.at(10, 1) = -0.5 * u * ( 1.0 + v ) * ( 1.0 + w );
213 dN.at(11, 1) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + w );
214 dN.at(12, 1) = -0.5 * u * ( 1.0 - v ) * ( 1.0 + w );
215 dN.at(13, 1) = -0.25 * ( 1.0 - v * v ) * ( 1.0 - w );
216 dN.at(14, 1) = -0.5 * u * ( 1.0 + v ) * ( 1.0 - w );
217 dN.at(15, 1) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - w );
218 dN.at(16, 1) = -0.5 * u * ( 1.0 - v ) * ( 1.0 - w );
219 dN.at(17, 1) = -0.25 * ( 1.0 - v ) * ( 1.0 - w * w );
220 dN.at(18, 1) = -0.25 * ( 1.0 + v ) * ( 1.0 - w * w );
221 dN.at(19, 1) = 0.25 * ( 1.0 + v ) * ( 1.0 - w * w );
222 dN.at(20, 1) = 0.25 * ( 1.0 - v ) * ( 1.0 - w * w );
223
224 dN.at(1, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 + w ) * ( 2.0 * v + u - w + 1.0 );
225 dN.at(2, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 + w ) * ( 2.0 * v - u + w - 1.0 );
226 dN.at(3, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 + w ) * ( 2.0 * v + u + w - 1.0 );
227 dN.at(4, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 + w ) * ( 2.0 * v - u - w + 1.0 );
228 dN.at(5, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 - w ) * ( 2.0 * v + u + w + 1.0 );
229 dN.at(6, 2) = 0.125 * ( 1.0 - u ) * ( 1.0 - w ) * ( 2.0 * v - u - w - 1.0 );
230 dN.at(7, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 - w ) * ( 2.0 * v + u - w - 1.0 );
231 dN.at(8, 2) = 0.125 * ( 1.0 + u ) * ( 1.0 - w ) * ( 2.0 * v - u + w + 1.0 );
232 dN.at(9, 2) = -0.5 * v * ( 1.0 - u ) * ( 1.0 + w );
233 dN.at(10, 2) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + w );
234 dN.at(11, 2) = -0.5 * v * ( 1.0 + u ) * ( 1.0 + w );
235 dN.at(12, 2) = -0.25 * ( 1.0 - u * u ) * ( 1.0 + w );
236 dN.at(13, 2) = -0.5 * v * ( 1.0 - u ) * ( 1.0 - w );
237 dN.at(14, 2) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - w );
238 dN.at(15, 2) = -0.5 * v * ( 1.0 + u ) * ( 1.0 - w );
239 dN.at(16, 2) = -0.25 * ( 1.0 - u * u ) * ( 1.0 - w );
240 dN.at(17, 2) = -0.25 * ( 1.0 - u ) * ( 1.0 - w * w );
241 dN.at(18, 2) = 0.25 * ( 1.0 - u ) * ( 1.0 - w * w );
242 dN.at(19, 2) = 0.25 * ( 1.0 + u ) * ( 1.0 - w * w );
243 dN.at(20, 2) = -0.25 * ( 1.0 + u ) * ( 1.0 - w * w );
244
245 dN.at(1, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 2.0 * w - u - v - 1.0 );
246 dN.at(2, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 2.0 * w - u + v - 1.0 );
247 dN.at(3, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 2.0 * w + u + v - 1.0 );
248 dN.at(4, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 2.0 * w + u - v - 1.0 );
249 dN.at(5, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 - v ) * ( 2.0 * w + u + v + 1.0 );
250 dN.at(6, 3) = 0.125 * ( 1.0 - u ) * ( 1.0 + v ) * ( 2.0 * w + u - v + 1.0 );
251 dN.at(7, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 + v ) * ( 2.0 * w - u - v + 1.0 );
252 dN.at(8, 3) = 0.125 * ( 1.0 + u ) * ( 1.0 - v ) * ( 2.0 * w - u + v + 1.0 );
253 dN.at(9, 3) = 0.25 * ( 1.0 - v * v ) * ( 1.0 - u );
254 dN.at(10, 3) = 0.25 * ( 1.0 - u * u ) * ( 1.0 + v );
255 dN.at(11, 3) = 0.25 * ( 1.0 - v * v ) * ( 1.0 + u );
256 dN.at(12, 3) = 0.25 * ( 1.0 - u * u ) * ( 1.0 - v );
257 dN.at(13, 3) = -0.25 * ( 1.0 - v * v ) * ( 1.0 - u );
258 dN.at(14, 3) = -0.25 * ( 1.0 - u * u ) * ( 1.0 + v );
259 dN.at(15, 3) = -0.25 * ( 1.0 - v * v ) * ( 1.0 + u );
260 dN.at(16, 3) = -0.25 * ( 1.0 - u * u ) * ( 1.0 - v );
261 dN.at(17, 3) = -0.5 * ( 1.0 - u ) * ( 1.0 - v ) * w;
262 dN.at(18, 3) = -0.5 * ( 1.0 - u ) * ( 1.0 + v ) * w;
263 dN.at(19, 3) = -0.5 * ( 1.0 + u ) * ( 1.0 + v ) * w;
264 dN.at(20, 3) = -0.5 * ( 1.0 + u ) * ( 1.0 - v ) * w;
265#endif
266}
267
268
269std::pair<double, FloatMatrixF<3,20>>
270FEI3dHexaQuad :: evaldNdx(const FloatArrayF<3> &lcoords, const FEICellGeometry &cellgeo)
271{
272 auto dNduvw = evaldNdxi(lcoords);
273 FloatMatrixF<3,20> coords;
274 for ( int i = 0; i < 20; i++ ) {
276 coords.setColumn(cellgeo.giveVertexCoordinates(i+1), i);
277 }
278 auto jacT = dotT(dNduvw, coords);
279 return {det(jacT), dot(inv(jacT), dNduvw)};
280}
281
282
283double
284FEI3dHexaQuad :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
285{
286#if 0
287 auto tmp = evaldNdx(lcoords, cellgeo);
288 answer = tmp.second;
289 return tmp.first;
290#else
291 FloatMatrix jacobianMatrix, inv, dNduvw, coords;
292 this->evaldNdxi(dNduvw, lcoords, cellgeo);
293 coords.resize( 3, dNduvw.giveNumberOfRows() );
294 for ( int i = 1; i <= dNduvw.giveNumberOfRows(); i++ ) {
295 coords.setColumn(cellgeo.giveVertexCoordinates(i), i);
296 }
297 jacobianMatrix.beProductOf(coords, dNduvw);
298 inv.beInverseOf(jacobianMatrix);
299
300 answer.beProductOf(dNduvw, inv);
301 return jacobianMatrix.giveDeterminant();
302#endif
303}
304
305void
306FEI3dHexaQuad :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
307{
308 FloatArray n;
309
310 this->evalN(n, lcoords, cellgeo);
311 answer.clear();
312 for ( int i = 1; i <= n.giveSize(); i++ ) {
313 answer.add( n.at(i), cellgeo.giveVertexCoordinates(i) );
314 }
315}
316
317double FEI3dHexaQuad :: giveCharacteristicLength(const FEICellGeometry &cellgeo) const
318{
319 const auto &n1 = cellgeo.giveVertexCoordinates(1);
320 const auto &n2 = cellgeo.giveVertexCoordinates(7);
322 return distance(n1, n2);
323}
324
325
326#define POINT_TOL 1.e-3
327
328int
329FEI3dHexaQuad :: global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo) const
330{
331 FloatArray res, delta, guess;
332 FloatMatrix jac;
333 double convergence_limit, error = 0.0;
334
335 // find a suitable convergence limit
336 convergence_limit = 1e-6 * this->giveCharacteristicLength(cellgeo);
337
338 // setup initial guess
339 answer.resize( gcoords.giveSize() );
340 answer.zero();
341
342 // apply Newton-Raphson to solve the problem
343 for ( int nite = 0; nite < 40; nite++ ) {
344 // compute the residual
345 this->local2global(guess, answer, cellgeo);
346 res.beDifferenceOf(gcoords, guess);
347
348 // check for convergence
349 error = res.computeNorm();
350 if ( error < convergence_limit ) {
351 break;
352 }
353
354 // compute the corrections
355 this->giveJacobianMatrixAt(jac, answer, cellgeo);
356 jac.solveForRhs(res, delta);
357
358 // update guess
359 answer.add(delta);
360 }
361 if ( error > convergence_limit ) { // Imperfect, could give false negatives.
362 //OOFEM_ERROR("no convergence after 10 iterations");
363 answer.zero();
364 return false;
365 }
366
367 // check limits for each local coordinate [-1,1] for quadrilaterals. (different for other elements, typically [0,1]).
368 bool inside = true;
369 for ( int i = 1; i <= answer.giveSize(); i++ ) {
370 if ( answer.at(i) < ( -1. - POINT_TOL ) ) {
371 answer.at(i) = -1.;
372 inside = false;
373 } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
374 answer.at(i) = 1.;
375 inside = false;
376 }
377 }
378
379 return inside;
380}
381
382void FEI3dHexaQuad :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
383{
384 double u = lcoords.at(1);
385 answer.resize(3);
386 answer.at(1) = 0.5 * ( u - 1. ) * u;
387 answer.at(2) = 0.5 * ( u + 1. ) * u;
388 answer.at(3) = 1. - u * u;
389}
390
391void FEI3dHexaQuad :: edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
392{
393 double u = lcoords.at(1);
394 const auto &eNodes = this->computeLocalEdgeMapping(iedge);
395 answer.clear();
396 answer.add( 0.5 * ( u - 1. ) * u, cellgeo.giveVertexCoordinates( eNodes.at(1) ) );
397 answer.add( 0.5 * ( u - 1. ) * u, cellgeo.giveVertexCoordinates( eNodes.at(2) ) );
398 answer.add( 1. - u * u, cellgeo.giveVertexCoordinates( eNodes.at(3) ) );
399}
400
401void FEI3dHexaQuad :: edgeEvaldNdx(FloatMatrix &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
402{
404 FloatArray dNdu;
405 double u = lcoords.at(1);
406 const auto &eNodes = this->computeLocalEdgeMapping(iedge);
407 dNdu.add( u - 0.5, cellgeo.giveVertexCoordinates( eNodes.at(1) ) );
408 dNdu.add( u + 0.5, cellgeo.giveVertexCoordinates( eNodes.at(2) ) );
409 dNdu.add( -2. * u, cellgeo.giveVertexCoordinates( eNodes.at(3) ) );
410 // Why matrix output?
411 answer.resize(3, 1);
412 answer.setColumn(dNdu, 1);
413}
414
415double FEI3dHexaQuad :: edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
416{
417 FloatArray dNdu;
418 double u = lcoords.at(1);
419 const auto &eNodes = this->computeLocalEdgeMapping(iedge);
420 dNdu.add( u - 0.5, cellgeo.giveVertexCoordinates( eNodes.at(1) ) );
421 dNdu.add( u + 0.5, cellgeo.giveVertexCoordinates( eNodes.at(2) ) );
422 dNdu.add( -2. * u, cellgeo.giveVertexCoordinates( eNodes.at(3) ) );
423 return dNdu.computeNorm();
424}
425
426IntArray
427FEI3dHexaQuad :: computeLocalEdgeMapping(int iedge) const
428{
429 if ( iedge == 1 ) {
430 return { 1, 2, 9};
431 } else if ( iedge == 2 ) {
432 return { 2, 3, 10};
433 } else if ( iedge == 3 ) {
434 return { 3, 4, 11};
435 } else if ( iedge == 4 ) {
436 return { 4, 1, 12};
437 } else if ( iedge == 5 ) {
438 return { 5, 6, 13};
439 } else if ( iedge == 6 ) {
440 return { 6, 7, 14};
441 } else if ( iedge == 7 ) {
442 return { 7, 8, 15};
443 } else if ( iedge == 8 ) {
444 return { 8, 5, 16};
445 } else if ( iedge == 9 ) {
446 return { 1, 5, 17};
447 } else if ( iedge == 10 ) {
448 return { 2, 6, 18};
449 } else if ( iedge == 11 ) {
450 return { 3, 7, 19};
451 } else if ( iedge == 12 ) {
452 return { 4, 8, 20};
453 } else {
454 throw std::range_error("invalid edge number");
455 }
456}
457
458void
459FEI3dHexaQuad :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
460{
461 double ksi = lcoords.at(1);
462 double eta = lcoords.at(2);
463
464 answer.resize(8);
465 answer.at(1) = ( 1. + ksi ) * ( 1. + eta ) * 0.25 * ( ksi + eta - 1. );
466 answer.at(2) = ( 1. - ksi ) * ( 1. + eta ) * 0.25 * ( -ksi + eta - 1. );
467 answer.at(3) = ( 1. - ksi ) * ( 1. - eta ) * 0.25 * ( -ksi - eta - 1. );
468 answer.at(4) = ( 1. + ksi ) * ( 1. - eta ) * 0.25 * ( ksi - eta - 1. );
469 answer.at(5) = 0.5 * ( 1. - ksi * ksi ) * ( 1. + eta );
470 answer.at(6) = 0.5 * ( 1. - ksi ) * ( 1. - eta * eta );
471 answer.at(7) = 0.5 * ( 1. - ksi * ksi ) * ( 1. - eta );
472 answer.at(8) = 0.5 * ( 1. + ksi ) * ( 1. - eta * eta );
473}
474
475void
476FEI3dHexaQuad :: surfaceEvaldNdx(FloatMatrix &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
477{
478 // Note, this must be in correct order, not just the correct nodes, therefore we must use snodes;
479 const auto &snodes = this->computeLocalSurfaceMapping(isurf);
480
481 FloatArray lcoords_hex;
482
484#if 1
485 if ( isurf == 1 ) { // surface 1 - nodes 1 4 3 2
486 lcoords_hex = Vec3(-lcoords.at(1), -lcoords.at(2), 1);
487 } else if ( isurf == 2 ) { // surface 2 - nodes 5 6 7 8
488 lcoords_hex = Vec3(-lcoords.at(2), -lcoords.at(1), -1);
489 } else if ( isurf == 3 ) { // surface 3 - nodes 1 2 6 5
490 lcoords_hex = Vec3(-1, -lcoords.at(1), lcoords.at(2));
491 } else if ( isurf == 4 ) { // surface 4 - nodes 2 3 7 6
492 lcoords_hex = Vec3(-lcoords.at(1), 1, lcoords.at(2));
493 } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7
494 lcoords_hex = Vec3(1, lcoords.at(1), lcoords.at(2));
495 } else if ( isurf == 6 ) { // surface 6 - nodes 4 1 5 8
496 lcoords_hex = Vec3(lcoords.at(1), -1, lcoords.at(2));
497 } else {
498 OOFEM_ERROR("wrong surface number (%d)", isurf);
499 }
500#else
502 if ( isurf == 1 ) { // surface 1 - nodes 3 4 8 7
503 lcoords_hex = Vec3(-1, lcoords.at(1), lcoords.at(2));
504 } else if ( isurf == 2 ) { // surface 2 - nodes 2 1 5 6
505 lcoords_hex = Vec3(1, lcoords.at(1), lcoords.at(2));
506 } else if ( isurf == 3 ) { // surface 3 - nodes 3 7 6 2
507 lcoords_hex = Vec3(lcoords.at(1), -1, lcoords.at(2));
508 } else if ( isurf == 4 ) { // surface 4 - nodes 4 8 5 1
509 lcoords_hex = Vec3(lcoords.at(1), 1, lcoords.at(2));
510 } else if ( isurf == 5 ) { // surface 5 - nodes 3 2 1 4
511 lcoords_hex = Vec3(lcoords.at(1), lcoords.at(2), -1);
512 } else if ( isurf == 6 ) { // surface 6 - nodes 7 6 5 8
513 lcoords_hex = Vec3(lcoords.at(1), lcoords.at(2), 1);
514 } else {
515 OOFEM_ERROR("wrong surface number (%d)", isurf);
516 }
517#endif
518
519 FloatMatrix fullB;
520 this->evaldNdx(fullB, lcoords_hex, cellgeo);
521 answer.resize(snodes.giveSize(), 3);
522 for ( int i = 1; i <= snodes.giveSize(); ++i ) {
523 for ( int j = 1; j <= 3; ++j ) {
524 answer.at(i, j) = fullB.at(snodes.at(i), j);
525 }
526 }
527}
528
529
530double
531FEI3dHexaQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
532{
533 FloatArray a, b, dNdksi(8), dNdeta(8);
534
535 const auto &snodes = this->computeLocalSurfaceMapping(isurf);
536
537 double ksi = lcoords.at(1);
538 double eta = lcoords.at(2);
539
540 dNdksi.at(1) = 0.25 * ( 1. + eta ) * ( 2.0 * ksi + eta );
541 dNdksi.at(2) = -0.25 * ( 1. + eta ) * ( -2.0 * ksi + eta );
542 dNdksi.at(3) = -0.25 * ( 1. - eta ) * ( -2.0 * ksi - eta );
543 dNdksi.at(4) = 0.25 * ( 1. - eta ) * ( 2.0 * ksi - eta );
544 dNdksi.at(5) = -ksi * ( 1. + eta );
545 dNdksi.at(6) = -0.5 * ( 1. - eta * eta );
546 dNdksi.at(7) = -ksi * ( 1. - eta );
547 dNdksi.at(8) = 0.5 * ( 1. - eta * eta );
548
549 dNdeta.at(1) = 0.25 * ( 1. + ksi ) * ( 2.0 * eta + ksi );
550 dNdeta.at(2) = 0.25 * ( 1. - ksi ) * ( 2.0 * eta - ksi );
551 dNdeta.at(3) = -0.25 * ( 1. - ksi ) * ( -2.0 * eta - ksi );
552 dNdeta.at(4) = -0.25 * ( 1. + ksi ) * ( -2.0 * eta + ksi );
553 dNdeta.at(5) = 0.5 * ( 1. - ksi * ksi );
554 dNdeta.at(6) = -eta * ( 1. - ksi );
555 dNdeta.at(7) = -0.5 * ( 1. - ksi * ksi );
556 dNdeta.at(8) = -eta * ( 1. + ksi );
557
558 for ( int i = 1; i <= 8; ++i ) {
559 a.add( dNdksi.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
560 b.add( dNdeta.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
561 }
562
563 answer.beVectorProductOf(a, b);
564 return answer.normalize_giveNorm();
565}
566
567void
568FEI3dHexaQuad :: surfaceLocal2global(FloatArray &answer, int isurf,
569 const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
570{
571 FloatArray n;
572
573 const auto &nodes = this->computeLocalSurfaceMapping(isurf);
574
575 this->surfaceEvalN(n, isurf, lcoords, cellgeo);
576
577 answer.clear();
578 for ( int i = 1; i <= n.giveSize(); ++i ) {
579 answer.add( n.at(i), cellgeo.giveVertexCoordinates( nodes.at(i) ) );
580 }
581}
582
583
584double
585FEI3dHexaQuad :: surfaceGiveTransformationJacobian(int isurf, const FloatArray &lcoords,
586 const FEICellGeometry &cellgeo) const
587{
588 FloatArray normal;
589 return this->surfaceEvalNormal(normal, isurf, lcoords, cellgeo);
590}
591
592
593IntArray
594FEI3dHexaQuad :: computeLocalSurfaceMapping(int isurf) const
595{
596 // the actual numbering has a positive normal pointing outwards from the element - (LSpace compatible)
597
598 if ( isurf == 1 ) {
599 return { 3, 2, 1, 4, 10, 9, 12, 11};
600 } else if ( isurf == 2 ) {
601 return { 7, 8, 5, 6, 15, 16, 13, 14};
602 } else if ( isurf == 3 ) {
603 return { 2, 6, 5, 1, 18, 13, 17, 9};
604 } else if ( isurf == 4 ) {
605 return { 3, 7, 6, 2, 19, 14, 18, 10};
606 } else if ( isurf == 5 ) {
607 return { 3, 4, 8, 7, 11, 20, 15, 19};
608 } else if ( isurf == 6 ) {
609 return { 4, 1, 5, 8, 12, 17, 16, 20};
610 } else {
611 throw std::range_error("invalid surface number");
612 }
613
614#if 0
615 // this commented numbering is symmetrical with respect to local coordinate axes
616
617 if ( isurf == 1 ) { // surface 1 - nodes 3 2 1 4 10 9 12 11
618 nodes.at(1) = 3;
619 nodes.at(2) = 2;
620 nodes.at(3) = 1;
621 nodes.at(4) = 4;
622 nodes.at(5) = 10;
623 nodes.at(6) = 9;
624 nodes.at(7) = 12;
625 nodes.at(8) = 11;
626 } else if ( iSurf == 2 ) { // surface 2 - nodes 7 6 5 8 14 13 16 15
627 nodes.at(1) = 7;
628 nodes.at(2) = 6;
629 nodes.at(3) = 5;
630 nodes.at(4) = 8;
631 nodes.at(5) = 14;
632 nodes.at(6) = 13;
633 nodes.at(7) = 16;
634 nodes.at(8) = 15;
635 } else if ( iSurf == 3 ) { // surface 3 - nodes 2 1 5 6 9 17 13 18
636 nodes.at(1) = 2;
637 nodes.at(2) = 1;
638 nodes.at(3) = 5;
639 nodes.at(4) = 6;
640 nodes.at(5) = 9;
641 nodes.at(6) = 17;
642 nodes.at(7) = 13;
643 nodes.at(8) = 18;
644 } else if ( isurf == 4 ) { // surface 4 - nodes 3 7 6 2 19 14 18 10
645 nodes.at(1) = 3;
646 nodes.at(2) = 7;
647 nodes.at(3) = 6;
648 nodes.at(4) = 2;
649 nodes.at(5) = 19;
650 nodes.at(6) = 14;
651 nodes.at(7) = 18;
652 nodes.at(8) = 10;
653 } else if ( isurf == 5 ) { // surface 5 - nodes 3 4 8 7 11 20 15 19
654 nodes.at(1) = 3;
655 nodes.at(2) = 4;
656 nodes.at(3) = 8;
657 nodes.at(4) = 7;
658 nodes.at(5) = 11;
659 nodes.at(6) = 20;
660 nodes.at(7) = 15;
661 nodes.at(8) = 19;
662 } else if ( iSurf == 6 ) { // surface 6 - nodes 4 8 5 1 20 16 17 12
663 nodes.at(1) = 4;
664 nodes.at(2) = 8;
665 nodes.at(3) = 5;
666 nodes.at(4) = 1;
667 nodes.at(5) = 20;
668 nodes.at(6) = 16;
669 nodes.at(7) = 17;
670 nodes.at(8) = 12;
671 } else {
672 OOFEM_ERROR("wrong surface number (%d)", isurf);
673 }
674#endif
675}
676
677
678void
679FEI3dHexaQuad :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
680// Returns the jacobian matrix J (x,y,z)/(ksi,eta,dzeta) of the receiver.
681{
682 FloatMatrix dNduvw, coords;
683 this->evaldNdxi(dNduvw, lcoords, cellgeo);
684 coords.resize( 3, dNduvw.giveNumberOfRows() );
685 for ( int i = 1; i <= dNduvw.giveNumberOfRows(); i++ ) {
686 coords.setColumn(cellgeo.giveVertexCoordinates(i), i);
687 }
688 jacobianMatrix.beProductOf(coords, dNduvw);
689}
690
691
692double
693FEI3dHexaQuad :: evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo) const
694{
695 const auto &fNodes = this->computeLocalSurfaceMapping(iEdge);
696
697 const auto &c1 = cellgeo.giveVertexCoordinates( fNodes.at(1) );
698 const auto &c2 = cellgeo.giveVertexCoordinates( fNodes.at(2) );
699 const auto &c3 = cellgeo.giveVertexCoordinates( fNodes.at(3) );
700 const auto &c4 = cellgeo.giveVertexCoordinates( fNodes.at(4) );
701 const auto &c5 = cellgeo.giveVertexCoordinates( fNodes.at(5) );
702 const auto &c6 = cellgeo.giveVertexCoordinates( fNodes.at(6) );
703 const auto &c7 = cellgeo.giveVertexCoordinates( fNodes.at(7) );
704 const auto &c8 = cellgeo.giveVertexCoordinates( fNodes.at(8) );
705
706 // Generated with Mathematica (rather unwieldy expression, tried to simplify it as good as possible, but it could probably be better)
707 return (
708 c1(2) * ( c2(1) * ( -3 * c3(0) - 3 * c4(0) - 12 * c5(0) + 14 * c6(0) + 14 * c8(0) ) +
709 c3(1) * ( 3 * c2(0) - 3 * c4(0) - 6 * c5(0) - 6 * c6(0) + 6 * c7(0) + 6 * c8(0) ) +
710 c4(1) * ( 3 * c2(0) + 3 * c3(0) - 14 * c5(0) - 14 * c7(0) + 12 * c8(0) ) +
711 c5(1) * ( 12 * c2(0) + 6 * c3(0) + 14 * c4(0) - 4 * c6(0) - 8 * c7(0) - 60 * c8(0) ) +
712 c6(1) * ( -14 * c2(0) + 6 * c3(0) + 4 * c5(0) + 12 * c7(0) - 8 * c8(0) ) +
713 c7(1) * ( -6 * c3(0) + 14 * c4(0) + 8 * c5(0) - 12 * c6(0) - 4 * c8(0) ) +
714 c8(1) * ( -14 * c2(0) - 6 * c3(0) - 12 * c4(0) + 60 * c5(0) + 8 * c6(0) + 4 * c7(0) ) ) +
715 c2(2) * ( c1(1) * ( 3 * c3(0) + 3 * c4(0) + 12 * c5(0) - 14 * c6(0) - 14 * c8(0) ) +
716 c3(1) * ( -3 * c1(0) - 3 * c4(0) + 14 * c5(0) - 12 * c6(0) + 14 * c7(0) ) +
717 c4(1) * ( -3 * c1(0) + 3 * c3(0) + 6 * c5(0) - 6 * c6(0) - 6 * c7(0) + 6 * c8(0) ) +
718 c5(1) * ( -12 * c1(0) - 14 * c3(0) - 6 * c4(0) + 60 * c6(0) + 8 * c7(0) + 4 * c8(0) ) +
719 c6(1) * ( 14 * c1(0) + 12 * c3(0) + 6 * c4(0) - 60 * c5(0) - 4 * c7(0) - 8 * c8(0) ) +
720 c7(1) * ( -14 * c3(0) + 6 * c4(0) - 8 * c5(0) + 4 * c6(0) + 12 * c8(0) ) +
721 c8(1) * ( 14 * c1(0) - 6 * c4(0) - 4 * c5(0) + 8 * c6(0) - 12 * c7(0) ) ) +
722 c3(2) * ( c1(1) * ( -3 * c2(0) + 3 * c4(0) + 6 * c5(0) + 6 * c6(0) - 6 * c7(0) - 6 * c8(0) ) +
723 c2(1) * ( 3 * c1(0) + 3 * c4(0) - 14 * c5(0) + 12 * c6(0) - 14 * c7(0) ) +
724 c4(1) * ( -3 * c1(0) - 3 * c2(0) + 14 * c6(0) - 12 * c7(0) + 14 * c8(0) ) +
725 c5(1) * ( -6 * c1(0) + 14 * c2(0) - 4 * c6(0) + 8 * c7(0) - 12 * c8(0) ) +
726 c6(1) * ( -6 * c1(0) - 12 * c2(0) - 14 * c4(0) + 4 * c5(0) + 60 * c7(0) + 8 * c8(0) ) +
727 c7(1) * ( 6 * c1(0) + 14 * c2(0) + 12 * c4(0) - 8 * c5(0) - 60 * c6(0) - 4 * c8(0) ) +
728 c8(1) * ( 6 * c1(0) - 14 * c4(0) + 12 * c5(0) - 8 * c6(0) + 4 * c7(0) ) ) +
729 c4(2) * ( c1(1) * ( -3 * c2(0) - 3 * c3(0) + 14 * c5(0) + 14 * c7(0) - 12 * c8(0) ) +
730 c2(1) * ( 3 * c1(0) - 3 * c3(0) - 6 * c5(0) + 6 * c6(0) + 6 * c7(0) - 6 * c8(0) ) +
731 c3(1) * ( 3 * c1(0) + 3 * c2(0) - 14 * c6(0) + 12 * c7(0) - 14 * c8(0) ) +
732 c5(1) * ( -14 * c1(0) + 6 * c2(0) + 12 * c6(0) - 8 * c7(0) + 4 * c8(0) ) +
733 c6(1) * ( -6 * c2(0) + 14 * c3(0) - 12 * c5(0) - 4 * c7(0) + 8 * c8(0) ) +
734 c7(1) * ( -14 * c1(0) - 6 * c2(0) - 12 * c3(0) + 8 * c5(0) + 4 * c6(0) + 60 * c8(0) ) +
735 c8(1) * ( 12 * c1(0) + 6 * c2(0) + 14 * c3(0) - 4 * c5(0) - 8 * c6(0) - 60 * c7(0) ) ) +
736 c5(2) * ( c1(1) * ( -12 * c2(0) - 6 * c3(0) - 14 * c4(0) + 4 * c6(0) + 8 * c7(0) + 60 * c8(0) ) +
737 c2(1) * ( 12 * c1(0) + 14 * c3(0) + 6 * c4(0) - 60 * c6(0) - 8 * c7(0) - 4 * c8(0) ) +
738 c3(1) * ( 6 * c1(0) - 14 * c2(0) + 4 * c6(0) - 8 * c7(0) + 12 * c8(0) ) +
739 c4(1) * ( 14 * c1(0) - 6 * c2(0) - 12 * c6(0) + 8 * c7(0) - 4 * c8(0) ) +
740 c6(1) * ( -4 * c1(0) + 60 * c2(0) - 4 * c3(0) + 12 * c4(0) - 32 * c7(0) - 32 * c8(0) ) +
741 c7(1) * ( -8 * c1(0) + 8 * c2(0) + 8 * c3(0) - 8 * c4(0) + 32 * c6(0) - 32 * c8(0) ) +
742 c8(1) * ( -60 * c1(0) + 4 * c2(0) - 12 * c3(0) + 4 * c4(0) + 32 * c6(0) + 32 * c7(0) ) ) +
743 c6(2) * ( c1(1) * ( 14 * c2(0) - 6 * c3(0) - 4 * c5(0) - 12 * c7(0) + 8 * c8(0) ) +
744 c2(1) * ( -14 * c1(0) - 12 * c3(0) - 6 * c4(0) + 60 * c5(0) + 4 * c7(0) + 8 * c8(0) ) +
745 c3(1) * ( 6 * c1(0) + 12 * c2(0) + 14 * c4(0) - 4 * c5(0) - 60 * c7(0) - 8 * c8(0) ) +
746 c4(1) * ( 6 * c2(0) - 14 * c3(0) + 12 * c5(0) + 4 * c7(0) - 8 * c8(0) ) +
747 c5(1) * ( 4 * c1(0) - 60 * c2(0) + 4 * c3(0) - 12 * c4(0) + 32 * c7(0) + 32 * c8(0) ) +
748 c7(1) * ( 12 * c1(0) - 4 * c2(0) + 60 * c3(0) - 4 * c4(0) - 32 * c5(0) - 32 * c8(0) ) +
749 c8(1) * ( -8 * c1(0) - 8 * c2(0) + 8 * c3(0) + 8 * c4(0) - 32 * c5(0) + 32 * c7(0) ) ) +
750 c7(2) * ( c1(1) * ( 6 * c3(0) - 14 * c4(0) - 8 * c5(0) + 12 * c6(0) + 4 * c8(0) ) +
751 c2(1) * ( 14 * c3(0) - 6 * c4(0) + 8 * c5(0) - 4 * c6(0) - 12 * c8(0) ) +
752 c3(1) * ( -6 * c1(0) - 14 * c2(0) - 12 * c4(0) + 8 * c5(0) + 60 * c6(0) + 4 * c8(0) ) +
753 c4(1) * ( 14 * c1(0) + 6 * c2(0) + 12 * c3(0) - 8 * c5(0) - 4 * c6(0) - 60 * c8(0) ) +
754 c5(1) * ( 8 * c1(0) - 8 * c2(0) - 8 * c3(0) + 8 * c4(0) - 32 * c6(0) + 32 * c8(0) ) +
755 c6(1) * ( -12 * c1(0) + 4 * c2(0) - 60 * c3(0) + 4 * c4(0) + 32 * c5(0) + 32 * c8(0) ) +
756 c8(1) * ( -4 * c1(0) + 12 * c2(0) - 4 * c3(0) + 60 * c4(0) - 32 * c5(0) - 32 * c6(0) ) ) +
757 c8(2) * ( c1(1) * ( 14 * c2(0) + 6 * c3(0) + 12 * c4(0) - 60 * c5(0) - 8 * c6(0) - 4 * c7(0) ) +
758 c2(1) * ( -14 * c1(0) + 6 * c4(0) + 4 * c5(0) - 8 * c6(0) + 12 * c7(0) ) +
759 c3(1) * ( -6 * c1(0) + 14 * c4(0) - 12 * c5(0) + 8 * c6(0) - 4 * c7(0) ) +
760 c4(1) * ( -12 * c1(0) - 6 * c2(0) - 14 * c3(0) + 4 * c5(0) + 8 * c6(0) + 60 * c7(0) ) +
761 c5(1) * ( 60 * c1(0) - 4 * c2(0) + 12 * c3(0) - 4 * c4(0) - 32 * c6(0) - 32 * c7(0) ) +
762 c6(1) * ( 8 * c1(0) + 8 * c2(0) - 8 * c3(0) - 8 * c4(0) + 32 * c5(0) - 32 * c7(0) ) +
763 c7(1) * ( 4 * c1(0) - 12 * c2(0) + 4 * c3(0) - 60 * c4(0) + 32 * c5(0) + 32 * c6(0) ) )
764 ) / 60.0;
765}
766
767
768std::unique_ptr<IntegrationRule>
769FEI3dHexaQuad :: giveIntegrationRule(int order, Element_Geometry_Type egt) const
770{
771 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
772 int points = iRule->getRequiredNumberOfIntegrationPoints(_Cube, order + 9);
773 iRule->SetUpPointsOnCube(points, _Unknown);
774 return std::move(iRule);
775}
776
777std::unique_ptr<IntegrationRule>
778FEI3dHexaQuad :: giveBoundaryIntegrationRule(int order, int boundary, Element_Geometry_Type egt) const
779{
780 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
781 int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 4);
782 iRule->SetUpPointsOnSquare(points, _Unknown);
783 return std::move(iRule);
784}
785} // end namespace oofem
oofem::FEI3dHexaQuad::computeLocalSurfaceMapping
IntArray computeLocalSurfaceMapping(int iSurf) const override
Definition fei3dhexaquad.C:594
oofem::FEI3dHexaQuad::computeLocalEdgeMapping
IntArray computeLocalEdgeMapping(int iedge) const override
Definition fei3dhexaquad.C:427
oofem::FEI3dHexaQuad::surfaceEvalN
void surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dhexaquad.C:459
oofem::FEI3dHexaQuad::giveCharacteristicLength
double giveCharacteristicLength(const FEICellGeometry &cellgeo) const
Definition fei3dhexaquad.C:317
oofem::FEI3dHexaQuad::local2global
void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dhexaquad.C:306
oofem::FEI3dHexaQuad::surfaceEvalNormal
double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dhexaquad.C:531
oofem::FEI3dHexaQuad::giveJacobianMatrixAt
void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dhexaquad.C:679
oofem::FEI3dHexaQuad::evaldNdxi
static FloatMatrixF< 3, 20 > evaldNdxi(const FloatArrayF< 3 > &lcoords)
Definition fei3dhexaquad.C:118
oofem::FEI3dHexaQuad::evalN
static FloatArrayF< 20 > evalN(const FloatArrayF< 3 > &lcoords)
Definition fei3dhexaquad.C:46
oofem::FEI3dHexaQuad::evaldNdx
static std::pair< double, FloatMatrixF< 3, 20 > > evaldNdx(const FloatArrayF< 3 > &lcoords, const FEICellGeometry &cellgeo)
Definition fei3dhexaquad.C:270
oofem::FEICellGeometry::giveVertexCoordinates
virtual const FloatArray giveVertexCoordinates(int i) const =0
oofem::FloatArray::computeNorm
double computeNorm() const
Definition floatarray.C:861
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::giveSize
Index giveSize() const
Returns the size of receiver.
Definition floatarray.h:261
oofem::FloatArray::beDifferenceOf
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Definition floatarray.C:403
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::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::giveNumberOfRows
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition floatmatrix.h:200
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_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::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