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