OOFEM 3.0
Loading...
Searching...
No Matches
fei3dwedgequad.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 "fei3dwedgequad.h"
36#include "floatarray.h"
37#include "floatmatrix.h"
38#include "floatarrayf.h"
39#include "floatmatrixf.h"
40#include "intarray.h"
41#include "gaussintegrationrule.h"
42
43namespace oofem {
44
45FloatArrayF<15>
46FEI3dWedgeQuad :: evalN(const FloatArrayF<3> &lcoords)
47{
48 double x = lcoords[0];
49 double y = lcoords[1];
50 double z = lcoords[2];
51 return {
52 0.5 * ( ( 1. - x - y ) * ( 2. * ( 1. - x - y ) - 1. ) * ( 1. - z ) - ( 1. - x - y ) * ( 1. - z * z ) ),
53 0.5 * ( x * ( 2. * x - 1. ) * ( 1. - z ) - x * ( 1. - z * z ) ),
54 0.5 * ( y * ( 2. * y - 1. ) * ( 1. - z ) - y * ( 1. - z * z ) ),
55 0.5 * ( ( 1. - x - y ) * ( 2. * ( 1. - x - y ) - 1. ) * ( 1. + z ) - ( 1. - x - y ) * ( 1. - z * z ) ),
56 0.5 * ( x * ( 2. * x - 1. ) * ( 1. + z ) - x * ( 1. - z * z ) ),
57 0.5 * ( y * ( 2. * y - 1. ) * ( 1. + z ) - y * ( 1. - z * z ) ),
58 2. * ( 1. - x - y ) * x * ( 1. - z ),
59 2. * x * y * ( 1. - z ),
60 2. * y * ( 1. - x - y ) * ( 1. - z ),
61 2. * x * ( 1. - x - y ) * ( 1. + z ),
62 2. * x * y * ( 1. + z ),
63 2. * y * ( 1. - x - y ) * ( 1. + z ),
64 ( 1. - x - y ) * ( 1. - z * z ),
65 x * ( 1. - z * z ),
66 y * ( 1. - z * z ),
67 };
68}
69
70void
71FEI3dWedgeQuad :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
72{
73 double x, y, z;
74 answer.resize(15);
75
76 x = lcoords.at(1);
77 y = lcoords.at(2);
78 z = lcoords.at(3);
79
80 answer.at(1) = 0.5 * ( ( 1. - x - y ) * ( 2. * ( 1. - x - y ) - 1. ) * ( 1. - z ) - ( 1. - x - y ) * ( 1. - z * z ) );
81 answer.at(2) = 0.5 * ( x * ( 2. * x - 1. ) * ( 1. - z ) - x * ( 1. - z * z ) );
82 answer.at(3) = 0.5 * ( y * ( 2. * y - 1. ) * ( 1. - z ) - y * ( 1. - z * z ) );
83 answer.at(4) = 0.5 * ( ( 1. - x - y ) * ( 2. * ( 1. - x - y ) - 1. ) * ( 1. + z ) - ( 1. - x - y ) * ( 1. - z * z ) );
84 answer.at(5) = 0.5 * ( x * ( 2. * x - 1. ) * ( 1. + z ) - x * ( 1. - z * z ) );
85 answer.at(6) = 0.5 * ( y * ( 2. * y - 1. ) * ( 1. + z ) - y * ( 1. - z * z ) );
86 answer.at(7) = 2. * ( 1. - x - y ) * x * ( 1. - z );
87 answer.at(8) = 2. * x * y * ( 1. - z );
88 answer.at(9) = 2. * y * ( 1. - x - y ) * ( 1. - z );
89 answer.at(10) = 2. * x * ( 1. - x - y ) * ( 1. + z );
90 answer.at(11) = 2. * x * y * ( 1. + z );
91 answer.at(12) = 2. * y * ( 1. - x - y ) * ( 1. + z );
92 answer.at(13) = ( 1. - x - y ) * ( 1. - z * z );
93 answer.at(14) = x * ( 1. - z * z );
94 answer.at(15) = y * ( 1. - z * z );
95}
96
97
98FloatMatrixF<3,15>
99FEI3dWedgeQuad :: evaldNdxi(const FloatArrayF<3> &lcoords)
100{
101 double x = lcoords[0];
102 double y = lcoords[1];
103 double z = lcoords[2];
104
105 return {
106 1. / 2. - ( z - 1. ) * ( x + y - 1. ) - z * z / 2. - ( ( z - 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2., // 1
107 1. / 2. - ( z - 1. ) * ( x + y - 1. ) - z * z / 2. - ( ( z - 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2.,
108 -z * ( x + y - 1. ) - ( ( 2. * x + 2. * y - 1. ) * ( x + y - 1. ) ) / 2.,
109 z * z / 2. - x * ( z - 1. ) - ( ( 2. * x - 1. ) * ( z - 1. ) ) / 2. - 1 / 2., // 2
110 0.,
111 x * z - ( x * ( 2. * x - 1. ) ) / 2.,
112 0., // 3
113 z * z / 2. - y * ( z - 1. ) - ( ( 2. * y - 1. ) * ( z - 1. ) ) / 2. - 1. / 2.,
114 y * z - ( y * ( 2 * y - 1. ) ) / 2.,
115 ( ( z + 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2. + ( z + 1. ) * ( x + y - 1. ) - z * z / 2. + 1. / 2., // 4
116 ( ( z + 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2. + ( z + 1. ) * ( x + y - 1. ) - z * z / 2. + 1. / 2.,
117 ( ( 2. * x + 2. * y - 1. ) * ( x + y - 1. ) ) / 2. - z * ( x + y - 1. ),
118 ( ( 2. * x - 1. ) * ( z + 1. ) ) / 2. + x * ( z + 1. ) + z * z / 2. - 1. / 2., 0., ( x * ( 2. * x - 1. ) ) / 2. + x * z, // 5
119 0., ( ( 2. * y - 1. ) * ( z + 1. ) ) / 2. + y * ( z + 1. ) + z * z / 2. - 1. / 2., ( y * ( 2. * y - 1. ) ) / 2. + y * z, // 6
120 ( z - 1. ) * ( 2. * x + 2. * y - 2. ) + 2. * x * ( z - 1. ), 2. * x * ( z - 1. ), x * ( 2. * x + 2. * y - 2. ), // 7
121 -2. * y * ( z - 1. ), -2. * x * ( z - 1. ), -2. * x * y, // 8
122 2. * y * ( z - 1. ), 2. * ( z - 1. ) * ( x + y - 1. ) + 2. * y * ( z - 1. ), 2. * y * ( x + y - 1. ), // 9
123 -2. * ( z + 1. ) * ( x + y - 1. ) - 2. * x * ( z + 1. ), -2. * x * ( z + 1. ), -2. * x * ( x + y - 1. ), // 10
124 2. * y * ( z + 1. ), 2. * x * ( z + 1. ), 2. * x * y, // 11
125 -2. * y * ( z + 1. ), -2. * ( z + 1. ) * ( x + y - 1. ) - 2. * y * ( z + 1. ), -2. * y * ( x + y - 1. ), // 12
126 z * z - 1., z * z - 1., 2. * z * ( x + y - 1. ), // 13
127 1. - z * z, 0., -2. * x * z, // 14
128 0., 1. - z * z, -2. * y * z, // 15
129 };
130}
131
132
133void
134FEI3dWedgeQuad :: evaldNdxi(FloatMatrix &dN, const FloatArray &lcoords, const FEICellGeometry &) const
135{
136 double x, y, z;
137 x = lcoords.at(1);
138 y = lcoords.at(2);
139 z = lcoords.at(3);
140
141 dN.resize(15, 3);
142
143 dN.at(1, 1) = 1. / 2. - ( z - 1. ) * ( x + y - 1. ) - z * z / 2. - ( ( z - 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2.;
144 dN.at(2, 1) = z * z / 2. - x * ( z - 1. ) - ( ( 2. * x - 1. ) * ( z - 1. ) ) / 2. - 1 / 2.;
145 dN.at(3, 1) = 0.;
146 dN.at(4, 1) = ( ( z + 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2. + ( z + 1. ) * ( x + y - 1. ) - z * z / 2. + 1. / 2.;
147 dN.at(5, 1) = ( ( 2. * x - 1. ) * ( z + 1. ) ) / 2. + x * ( z + 1. ) + z * z / 2. - 1. / 2.;
148 dN.at(6, 1) = 0.;
149 dN.at(7, 1) = ( z - 1. ) * ( 2. * x + 2. * y - 2. ) + 2. * x * ( z - 1. );
150 dN.at(8, 1) = -2. * y * ( z - 1. );
151 dN.at(9, 1) = 2. * y * ( z - 1. );
152 dN.at(10, 1) = -2. * ( z + 1. ) * ( x + y - 1. ) - 2. * x * ( z + 1. );
153 dN.at(11, 1) = 2. * y * ( z + 1. );
154 dN.at(12, 1) = -2. * y * ( z + 1. );
155 dN.at(13, 1) = z * z - 1.;
156 dN.at(14, 1) = 1. - z * z;
157 dN.at(15, 1) = 0.;
158
159 dN.at(1, 2) = 1. / 2. - ( z - 1. ) * ( x + y - 1. ) - z * z / 2. - ( ( z - 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2.;
160 dN.at(2, 2) = 0.;
161 dN.at(3, 2) = z * z / 2. - y * ( z - 1. ) - ( ( 2. * y - 1. ) * ( z - 1. ) ) / 2. - 1. / 2.;
162 dN.at(4, 2) = ( ( z + 1. ) * ( 2. * x + 2. * y - 1. ) ) / 2. + ( z + 1. ) * ( x + y - 1. ) - z * z / 2. + 1. / 2.;
163 dN.at(5, 2) = 0.;
164 dN.at(6, 2) = ( ( 2. * y - 1. ) * ( z + 1. ) ) / 2. + y * ( z + 1. ) + z * z / 2. - 1. / 2.;
165 dN.at(7, 2) = 2. * x * ( z - 1. );
166 dN.at(8, 2) = -2. * x * ( z - 1. );
167 dN.at(9, 2) = 2. * ( z - 1. ) * ( x + y - 1. ) + 2. * y * ( z - 1. );
168 dN.at(10, 2) = -2. * x * ( z + 1. );
169 dN.at(11, 2) = 2. * x * ( z + 1. );
170 dN.at(12, 2) = -2. * ( z + 1. ) * ( x + y - 1. ) - 2. * y * ( z + 1. );
171 dN.at(13, 2) = z * z - 1.;
172 dN.at(14, 2) = 0.;
173 dN.at(15, 2) = 1. - z * z;
174
175 dN.at(1, 3) = -z * ( x + y - 1. ) - ( ( 2. * x + 2. * y - 1. ) * ( x + y - 1. ) ) / 2.;
176 dN.at(2, 3) = x * z - ( x * ( 2. * x - 1. ) ) / 2.;
177 dN.at(3, 3) = y * z - ( y * ( 2 * y - 1. ) ) / 2.;
178 dN.at(4, 3) = ( ( 2. * x + 2. * y - 1. ) * ( x + y - 1. ) ) / 2. - z * ( x + y - 1. );
179 dN.at(5, 3) = ( x * ( 2. * x - 1. ) ) / 2. + x * z;
180 dN.at(6, 3) = ( y * ( 2. * y - 1. ) ) / 2. + y * z;
181 dN.at(7, 3) = x * ( 2. * x + 2. * y - 2. );
182 dN.at(8, 3) = -2. * x * y;
183 dN.at(9, 3) = 2. * y * ( x + y - 1. );
184 dN.at(10, 3) = -2. * x * ( x + y - 1. );
185 dN.at(11, 3) = 2. * x * y;
186 dN.at(12, 3) = -2. * y * ( x + y - 1. );
187 dN.at(13, 3) = 2. * z * ( x + y - 1. );
188 dN.at(14, 3) = -2. * x * z;
189 dN.at(15, 3) = -2. * y * z;
190}
191
192
193std::pair<double, FloatMatrixF<3,15>>
194FEI3dWedgeQuad :: evaldNdx(const FloatArrayF<3> &lcoords, const FEICellGeometry &cellgeo)
195{
196 auto dNduvw = evaldNdxi(lcoords);
197 FloatMatrixF<3,15> coords;
198 for ( int i = 0; i < 15; i++ ) {
199 coords.setColumn(cellgeo.giveVertexCoordinates(i+1), i);
200 }
201 auto jacT = dotT(dNduvw, coords);
202 return {det(jacT), dot(inv(jacT), dNduvw)};
203}
204
205
206double
207FEI3dWedgeQuad :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
208{
209 FloatMatrix jacobianMatrix, inv, dNduvw, coords;
210 this->evaldNdxi(dNduvw, lcoords, cellgeo);
211 coords.resize(3, 15);
212 for ( int i = 1; i <= 15; i++ ) {
213 coords.setColumn(cellgeo.giveVertexCoordinates(i), i);
214 }
215 jacobianMatrix.beProductOf(coords, dNduvw);
216 inv.beInverseOf(jacobianMatrix);
217
218 answer.beProductOf(dNduvw, inv);
219 return jacobianMatrix.giveDeterminant();
220}
221
222
223void
224FEI3dWedgeQuad :: giveLocalNodeCoords(FloatMatrix &answer, const Element_Geometry_Type) const
225{
226
227 answer.resize(3,15);
228 answer.at(1,1) = 1.0; answer.at(2,1) = 0.0; answer.at(3,1) = -1.0;
229 answer.at(1,2) = 0.0; answer.at(2,2) = 1.0; answer.at(3,2) = -1.0;
230 answer.at(1,3) = 0.0; answer.at(2,3) = 0.0; answer.at(3,3) = -1.0;
231 answer.at(1,4) = 1.0; answer.at(2,4) = 0.0; answer.at(3,4) = 1.0;
232 answer.at(1,5) = 0.0; answer.at(2,5) = 1.0; answer.at(3,5) = 1.0;
233 answer.at(1,6) = 0.0; answer.at(2,6) = 0.0; answer.at(3,6) = 1.0;
234 answer.at(1,7) = 0.5; answer.at(2,7) = 0.5; answer.at(3,7) = -1.0;
235 answer.at(1,8) = 0.0; answer.at(2,8) = 0.5; answer.at(3,8) = -1.0;
236 answer.at(1,9) = 0.5; answer.at(2,9) = 0.0; answer.at(3,9) = -1.0;
237 answer.at(1,10) = 0.5; answer.at(2,10) = 0.5; answer.at(3,10) = 1.0;
238 answer.at(1,11) = 0.0; answer.at(2,11) = 0.5; answer.at(3,11) = 1.0;
239 answer.at(1,12) = 0.5; answer.at(2,12) = 0.0; answer.at(3,12) = 1.0;
240 answer.at(1,13) = 1.0; answer.at(2,13) = 0.0; answer.at(3,13) = 0.0;
241 answer.at(1,14) = 0.0; answer.at(2,14) = 1.0; answer.at(3,14) = 0.0;
242 answer.at(1,15) = 0.0; answer.at(2,15) = 0.0; answer.at(3,15) = 0.0;
243
244}
245
246
247void
248FEI3dWedgeQuad :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
249{
250 FloatArray n;
251 this->evalN(n, lcoords, cellgeo);
252
253 answer.resize(3);
254 answer.zero();
255 for ( int i = 1; i <= 15; i++ ) {
256 answer.at(1) += n.at(i) * cellgeo.giveVertexCoordinates(i).at(1);
257 answer.at(2) += n.at(i) * cellgeo.giveVertexCoordinates(i).at(2);
258 answer.at(3) += n.at(i) * cellgeo.giveVertexCoordinates(i).at(3);
259 }
260}
261
262#define POINT_TOL 1.e-3
263
264int
265FEI3dWedgeQuad :: global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo) const
266{
267 FloatArray res, delta, guess;
268 FloatMatrix jac;
269 double convergence_limit, error = 0.0;
270
271 // find a suitable convergence limit
272 convergence_limit = 1e-6 * this->giveCharacteristicLength(cellgeo);
273
274 // setup initial guess
275 answer.resize( gcoords.giveSize() );
276 answer.zero();
277
278 // apply Newton-Raphson to solve the problem
279 for ( int nite = 0; nite < 40; nite++ ) {
280 // compute the residual
281 this->local2global(guess, answer, cellgeo);
282 res.beDifferenceOf(gcoords, guess);
283
284 // check for convergence
285 error = res.computeNorm();
286 if ( error < convergence_limit ) {
287 break;
288 }
289
290 // compute the corrections
291 this->giveJacobianMatrixAt(jac, answer, cellgeo);
292 jac.solveForRhs(res, delta);
293
294 // update guess
295 answer.add(delta);
296 }
297 if ( error > convergence_limit ) { // Imperfect, could give false negatives.
298 //OOFEM_ERROR("no convergence after 10 iterations");
299 answer.zero();
300 return false;
301 }
302
303 // check limits for each local coordinate [-1,1] for quadrilaterals. (different for other elements, typically [0,1]).
304 bool inside = true;
305 for ( int i = 1; i <= answer.giveSize(); i++ ) {
306 if ( answer.at(i) < ( -1. - POINT_TOL ) ) {
307 answer.at(i) = -1.;
308 inside = false;
309 } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
310 answer.at(i) = 1.;
311 inside = false;
312 }
313 }
314
315 return inside;
316}
317
318double FEI3dWedgeQuad :: giveCharacteristicLength(const FEICellGeometry &cellgeo) const
319{
320 const auto &n1 = cellgeo.giveVertexCoordinates(1);
321 const auto &n2 = cellgeo.giveVertexCoordinates(4);
323 return distance(n1, n2);
324}
325
326
327
328double
329FEI3dWedgeQuad :: giveTransformationJacobian(const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
330{
331 FloatMatrix jacobianMatrix;
332
333 this->giveJacobianMatrixAt(jacobianMatrix, lcoords, cellgeo);
334 return jacobianMatrix.giveDeterminant();
335}
336
337
338void
339FEI3dWedgeQuad :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
340// Returns the jacobian matrix J (x,y,z)/(ksi,eta,dzeta) of the receiver.
341{
342 FloatMatrix dNduvw, coords;
343 this->evaldNdxi(dNduvw, lcoords, cellgeo);
344 coords.resize(3, 15);
345 for ( int i = 1; i <= 15; i++ ) {
346 coords.setColumn(cellgeo.giveVertexCoordinates(i), i);
347 }
348 jacobianMatrix.beProductOf(coords, dNduvw);
349}
350
351
352void FEI3dWedgeQuad :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
353{
354 double ksi = lcoords.at(1);
355 answer.resize(3);
356 answer.at(1) = ksi * ( ksi - 1. ) * 0.5;
357 answer.at(2) = ksi * ( 1. + ksi ) * 0.5;
358 answer.at(3) = ( 1. - ksi * ksi );
359}
360
361
362void FEI3dWedgeQuad :: edgeEvaldNdx(FloatMatrix &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
363{
364 OOFEM_ERROR("not implemented");
365}
366
367
368void FEI3dWedgeQuad :: edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
369{
370 FloatArray n;
371
372 const auto &nodes = this->computeLocalEdgeMapping(iedge);
373 this->edgeEvalN(n, iedge, lcoords, cellgeo);
374
375 answer.clear();
376 for ( int i = 1; i <= n.giveSize(); ++i ) {
377 answer.add( n.at(i), cellgeo.giveVertexCoordinates( nodes.at(i) ) );
378 }
379}
380
381
382IntArray
383FEI3dWedgeQuad :: computeLocalEdgeMapping(int iedge) const
384{
385 if ( iedge == 1 ) {
386 return {1, 2, 7};
387 } else if ( iedge == 2 ) {
388 return {2, 3, 8};
389 } else if ( iedge == 3 ) {
390 return {3, 1, 9};
391 } else if ( iedge == 4 ) {
392 return {4, 5, 10};
393 } else if ( iedge == 5 ) {
394 return {5, 6, 11};
395 } else if ( iedge == 6 ) {
396 return {6, 4, 12};
397 } else if ( iedge == 7 ) {
398 return {1, 4, 13};
399 } else if ( iedge == 8 ) {
400 return {2, 5, 14};
401 } else if ( iedge == 9 ) {
402 return {3, 6, 15};
403 } else {
404 throw std::range_error("invalid edge number");
405 //return {};
406 }
407}
408
409
410double FEI3dWedgeQuad :: edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
411{
412 OOFEM_ERROR("not implemented");
413 // return 0.0;
414}
415
416
417void
418FEI3dWedgeQuad :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
419{
420 if ( isurf <= 2 ) {
421 double l1, l2, l3;
422 l1 = lcoords.at(1);
423 l2 = lcoords.at(2);
424 l3 = 1.0 - l1 - l2;
425
426 answer.resize(6);
427 answer.at(1) = ( 2. * l1 - 1. ) * l1;
428 answer.at(2) = ( 2. * l2 - 1. ) * l2;
429 answer.at(3) = ( 2. * l3 - 1. ) * l3;
430 answer.at(4) = 4. * l1 * l2;
431 answer.at(5) = 4. * l2 * l3;
432 answer.at(6) = 4. * l3 * l1;
433 } else {
434 double ksi = lcoords.at(1);
435 double eta = lcoords.at(2);
436
437 answer.resize(8);
438 answer.at(1) = ( 1. + ksi ) * ( 1. + eta ) * 0.25 * ( ksi + eta - 1. );
439 answer.at(2) = ( 1. - ksi ) * ( 1. + eta ) * 0.25 * ( -ksi + eta - 1. );
440 answer.at(3) = ( 1. - ksi ) * ( 1. - eta ) * 0.25 * ( -ksi - eta - 1. );
441 answer.at(4) = ( 1. + ksi ) * ( 1. - eta ) * 0.25 * ( ksi - eta - 1. );
442 answer.at(5) = 0.5 * ( 1. - ksi * ksi ) * ( 1. + eta );
443 answer.at(6) = 0.5 * ( 1. - ksi ) * ( 1. - eta * eta );
444 answer.at(7) = 0.5 * ( 1. - ksi * ksi ) * ( 1. - eta );
445 answer.at(8) = 0.5 * ( 1. + ksi ) * ( 1. - eta * eta );
446 }
447}
448
449
450void
451FEI3dWedgeQuad :: surfaceLocal2global(FloatArray &answer, int isurf,
452 const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
453{
454 FloatArray n;
455
456 const auto &nodes = this->computeLocalSurfaceMapping(isurf);
457 this->surfaceEvalN(n, isurf, lcoords, cellgeo);
458
459 answer.clear();
460 for ( int i = 1; i <= n.giveSize(); ++i ) {
461 answer.add( n.at(i), cellgeo.giveVertexCoordinates( nodes.at(i) ) );
462 }
463}
464
465
466IntArray
467FEI3dWedgeQuad :: computeLocalSurfaceMapping(int isurf) const
468{
469 if ( isurf == 1 ) {
470 return {1, 2, 3, 7, 8, 9};
471 } else if ( isurf == 2 ) {
472 return {4, 5, 6, 10, 11, 12};
473 } else if ( isurf == 3 ) {
474 return {1, 2, 5, 4, 7, 14, 10, 13};
475 } else if ( isurf == 4 ) {
476 return {2, 3, 6, 5, 8, 15, 11, 14};
477 } else if ( isurf == 5 ) {
478 return {3, 1, 4, 6, 9, 13, 12, 15};
479 } else {
480 OOFEM_ERROR("Surface %d doesn't exist.\n", isurf);
481 // return {};
482 }
483}
484
485
486double
487FEI3dWedgeQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
488{
489 FloatArray a, b, dNdksi, dNdeta;
490
491 const auto &snodes = this->computeLocalSurfaceMapping(isurf);
492
493 if ( snodes.giveSize() == 6 ) {
494 double l1, l2, l3;
495 l1 = lcoords.at(1);
496 l2 = lcoords.at(2);
497 l3 = 1.0 - l1 - l2;
498
499 dNdksi.resize(6);
500 dNdksi(0) = 4.0 * l1 - 1.0;
501 dNdksi(1) = 0.0;
502 dNdksi(2) = -1.0 * ( 4.0 * l3 - 1.0 );
503 dNdksi(3) = 4.0 * l2;
504 dNdksi(4) = -4.0 * l2;
505 dNdksi(5) = 4.0 * l3 - 4.0 * l1;
506
507 dNdeta.resize(6);
508 dNdeta(0) = 0.0;
509 dNdeta(1) = 4.0 * l2 - 1.0;
510 dNdeta(2) = -1.0 * ( 4.0 * l3 - 1.0 );
511 dNdeta(3) = 4.0 * l1;
512 dNdeta(4) = 4.0 * l3 - 4.0 * l2;
513 dNdeta(5) = -4.0 * l1;
514 } else {
515 double ksi, eta;
516 ksi = lcoords.at(1);
517 eta = lcoords.at(2);
518
519 dNdksi.resize(8);
520 dNdksi.at(1) = 0.25 * ( 1. + eta ) * ( 2.0 * ksi + eta );
521 dNdksi.at(2) = -0.25 * ( 1. + eta ) * ( -2.0 * ksi + eta );
522 dNdksi.at(3) = -0.25 * ( 1. - eta ) * ( -2.0 * ksi - eta );
523 dNdksi.at(4) = 0.25 * ( 1. - eta ) * ( 2.0 * ksi - eta );
524 dNdksi.at(5) = -ksi * ( 1. + eta );
525 dNdksi.at(6) = -0.5 * ( 1. - eta * eta );
526 dNdksi.at(7) = -ksi * ( 1. - eta );
527 dNdksi.at(8) = 0.5 * ( 1. - eta * eta );
528
529 dNdeta.resize(8);
530 dNdeta.at(1) = 0.25 * ( 1. + ksi ) * ( 2.0 * eta + ksi );
531 dNdeta.at(2) = 0.25 * ( 1. - ksi ) * ( 2.0 * eta - ksi );
532 dNdeta.at(3) = -0.25 * ( 1. - ksi ) * ( -2.0 * eta - ksi );
533 dNdeta.at(4) = -0.25 * ( 1. + ksi ) * ( -2.0 * eta + ksi );
534 dNdeta.at(5) = 0.5 * ( 1. - ksi * ksi );
535 dNdeta.at(6) = -eta * ( 1. - ksi );
536 dNdeta.at(7) = -0.5 * ( 1. - ksi * ksi );
537 dNdeta.at(8) = -eta * ( 1. + ksi );
538 }
539
540 for ( int i = 1; i <= snodes.giveSize(); ++i ) {
541 a.add( dNdksi.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
542 b.add( dNdeta.at(i), cellgeo.giveVertexCoordinates( snodes.at(i) ) );
543 }
544
545 answer.beVectorProductOf(a, b);
546 return answer.normalize_giveNorm();
547}
548
549
550double
551FEI3dWedgeQuad :: surfaceGiveTransformationJacobian(int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
552{
553 FloatArray normal;
554 return this->surfaceEvalNormal(normal, isurf, lcoords, cellgeo);
555}
556
557
558std::unique_ptr<IntegrationRule>
559FEI3dWedgeQuad :: giveIntegrationRule(int order, Element_Geometry_Type egt) const
560{
561 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
563 //int points = iRule->getRequiredNumberOfIntegrationPoints(_Wedge, order);
564 OOFEM_WARNING("Warning.. ignoring 'order' argument: FIXME");
565 int pointsZeta = 1;
566 int pointsTriangle = 1;
567 iRule->SetUpPointsOnWedge(pointsTriangle, pointsZeta, _Unknown);
568 return std::move(iRule);
569}
570
571std::unique_ptr<IntegrationRule>
572FEI3dWedgeQuad :: giveBoundaryIntegrationRule(int order, int boundary, Element_Geometry_Type egt) const
573{
574 auto iRule = std::make_unique<GaussIntegrationRule>(1, nullptr);
575 if ( boundary <= 2 ) {
576 int points = iRule->getRequiredNumberOfIntegrationPoints(_Triangle, order + 2);
577 iRule->SetUpPointsOnTriangle(points, _Unknown);
578 } else {
580 int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 2);
581 iRule->SetUpPointsOnSquare(points, _Unknown);
582 }
583 return std::move(iRule);
584}
585} // end namespace oofem
oofem::FEI3dWedgeQuad::surfaceEvalNormal
double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dwedgequad.C:487
oofem::FEI3dWedgeQuad::surfaceEvalN
void surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dwedgequad.C:418
oofem::FEI3dWedgeQuad::evalN
static FloatArrayF< 15 > evalN(const FloatArrayF< 3 > &lcoords)
Definition fei3dwedgequad.C:46
oofem::FEI3dWedgeQuad::local2global
void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dwedgequad.C:248
oofem::FEI3dWedgeQuad::evaldNdxi
static FloatMatrixF< 3, 15 > evaldNdxi(const FloatArrayF< 3 > &lcoords)
Definition fei3dwedgequad.C:99
oofem::FEI3dWedgeQuad::computeLocalEdgeMapping
IntArray computeLocalEdgeMapping(int iedge) const override
Definition fei3dwedgequad.C:383
oofem::FEI3dWedgeQuad::computeLocalSurfaceMapping
IntArray computeLocalSurfaceMapping(int iSurf) const override
Definition fei3dwedgequad.C:467
oofem::FEI3dWedgeQuad::giveJacobianMatrixAt
void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dwedgequad.C:339
oofem::FEI3dWedgeQuad::edgeEvalN
void edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const override
Definition fei3dwedgequad.C:352
oofem::FEI3dWedgeQuad::giveCharacteristicLength
double giveCharacteristicLength(const FEICellGeometry &cellgeo) const
Definition fei3dwedgequad.C:318
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::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_WARNING
#define OOFEM_WARNING(...)
Definition error.h:80
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
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