OOFEM  2.4 OOFEM.org - Object Oriented Finite Element Solver
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 - 2013 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
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
36 #include "mathfem.h"
37 #include "floatmatrix.h"
38 #include "floatarray.h"
39 #include "gaussintegrationrule.h"
40
41 namespace oofem {
42 double
44 {
45  // Use linear approximation.
46
47  double x1, x2, x3, x4; //, x5, x6, x7, x8, x9, x10;
48  double y1, y2, y3, y4; //, y5, y6, y7, y8, y9, y10;
49  double z1, z2, z3, z4; //, z5, z6, z7, z8, z9, z10;
50
51  x1 = cellgeo.giveVertexCoordinates(1)->at(1);
52  y1 = cellgeo.giveVertexCoordinates(1)->at(2);
53  z1 = cellgeo.giveVertexCoordinates(1)->at(3);
54  x2 = cellgeo.giveVertexCoordinates(2)->at(1);
55  y2 = cellgeo.giveVertexCoordinates(2)->at(2);
56  z2 = cellgeo.giveVertexCoordinates(2)->at(3);
57  x3 = cellgeo.giveVertexCoordinates(3)->at(1);
58  y3 = cellgeo.giveVertexCoordinates(3)->at(2);
59  z3 = cellgeo.giveVertexCoordinates(3)->at(3);
60  x4 = cellgeo.giveVertexCoordinates(4)->at(1);
61  y4 = cellgeo.giveVertexCoordinates(4)->at(2);
62  z4 = cellgeo.giveVertexCoordinates(4)->at(3);
63 #if 0
64  x5=cellgeo.giveVertexCoordinates(5)->at(1); y5=cellgeo.giveVertexCoordinates(5)->at(2); z5=cellgeo.giveVertexCoordinates(5)->at(3);
65  x6=cellgeo.giveVertexCoordinates(6)->at(1); y6=cellgeo.giveVertexCoordinates(6)->at(2); z6=cellgeo.giveVertexCoordinates(6)->at(3);
66  x7=cellgeo.giveVertexCoordinates(7)->at(1); y7=cellgeo.giveVertexCoordinates(7)->at(2); z7=cellgeo.giveVertexCoordinates(7)->at(3);
67  x8=cellgeo.giveVertexCoordinates(8)->at(1); y8=cellgeo.giveVertexCoordinates(8)->at(2); z8=cellgeo.giveVertexCoordinates(8)->at(3);
68  x9=cellgeo.giveVertexCoordinates(9)->at(1); y9=cellgeo.giveVertexCoordinates(9)->at(2); z9=cellgeo.giveVertexCoordinates(9)->at(3);
69  x10=cellgeo.giveVertexCoordinates(10)->at(1); y10=cellgeo.giveVertexCoordinates(10)->at(2); z10=cellgeo.giveVertexCoordinates(10)->at(3);
70
71  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
72  + 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
73  + 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
74  + 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
75  + 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
76  + - 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
77  + - 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
78  + - 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
79  + - 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
80  + 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
81  + - 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
82  + - 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
83  + - 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
84  + 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
85  + 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
86  + 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
87  + 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
88  + - 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 -
89  + 2*x10*y4*z1 + 4*x1*y6*z9 - 4*x1*y7*z8 + 4*x1*y8*z7 - 4*x1*y9*z6
90  + - 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
91  + 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
92  + 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
93  + - 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
94  + - 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
95  + 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
96  + 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
97  + 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
98  + - 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
99  + - 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
100  + 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
101  + - 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
102  + 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
103  + - 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
104  + - 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
105  + 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
106  + - 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
107  + - 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
108  + 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
109  + - 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
110  + 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
111  + - 16*x8*y10*z9 - 16*x9*y8*z10 + 16*x9*y10*z8 + 16*x10*y8*z9 - 16*x10*y9*z8;
112
113  //printf("Q Area=%f\n", area);
114 #endif
115  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 -
116  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;
117
118  //printf("L Area=%f\n", area);
119
120  area = area / 6.0;
121  return area;
122 }
123
124 void
126 {
127  double x1 = lcoords(0);
128  double x2 = lcoords(1);
129  double x3 = lcoords(2);
130  double x4 = 1.0 - x1 - x2 - x3;
131
133  answer(0) = x1 * ( 2 * x1 - 1 );
134  answer(1) = x2 * ( 2 * x2 - 1 );
135  answer(2) = x3 * ( 2 * x3 - 1 );
136  answer(3) = x4 * ( 2 * x4 - 1 );
137
138  answer(4) = 4 * x1 * x2;
139  answer(5) = 4 * x2 * x3;
140  answer(6) = 4 * x3 * x1;
141  answer(7) = 4 * x1 * x4;
142  answer(8) = 4 * x2 * x4;
143  answer(9) = 4 * x3 * x4;
144 }
145
146 double
148 {
149  FloatMatrix jacobianMatrix, inv, dNduvw, coords;
150  this->evaldNdxi(dNduvw, lcoords, cellgeo);
151  coords.resize( 3, dNduvw.giveNumberOfRows() );
152  for ( int i = 1; i <= dNduvw.giveNumberOfRows(); i++ ) {
153  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
154  }
155  jacobianMatrix.beProductOf(coords, dNduvw);
156  inv.beInverseOf(jacobianMatrix);
157
159  return jacobianMatrix.giveDeterminant();
160 }
161
162 void
163 FEI3dTetQuad :: evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords , const FEICellGeometry &cellgeo)
164 {
165  double x1 = lcoords(0);
166  double x2 = lcoords(1);
167  double x3 = lcoords(2);
168  double x4 = 1.0 - x1 - x2 - x3;
169
171
172  // dNj/dx1
173  answer(0, 0) = 4 * x1 - 1;
176  answer(3, 0) = -4 * x4 + 1;
177  answer(4, 0) = 4 * x2;
179  answer(6, 0) = 4 * x3;
180  answer(7, 0) = 4 * ( x4 - x1 );
181  answer(8, 0) = -4 * x2;
182  answer(9, 0) = -4 * x3;
183
184  // dNj/dx2
186  answer(1, 1) = 4 * x2 - 1;
188  answer(3, 1) = -4 * x4 + 1;
189  answer(4, 1) = 4 * x1;
190  answer(5, 1) = 4 * x3;
192  answer(7, 1) = -4 * x1;
193  answer(8, 1) = 4 * ( x4 - x2 );
194  answer(9, 1) = -4 * x3;
195
196  // dNj/dx3
199  answer(2, 2) = 4 * x3 - 1;
200  answer(3, 2) = -4 * x4 + 1;
202  answer(5, 2) = 4 * x2;
203  answer(6, 2) = 4 * x1;
204  answer(7, 2) = -4 * x1;
205  answer(8, 2) = -4 * x2;
206  answer(9, 2) = 4 * ( x4 - x3 );
207 }
208
209
210 void
212 {
213  FloatArray N;
214  this->evalN(N, lcoords, cellgeo);
216  for ( int i = 1; i <= N.giveSize(); i++ ) {
218  }
219 }
220
221 #define POINT_TOL 1e-6
222
223 int
225 {
226  FloatArray res, delta, guess, lcoords_guess;
227  FloatMatrix jac;
228  double convergence_limit, error = 0.0;
229
230  // find a suitable convergence limit
231  convergence_limit = 1e-6 * this->giveCharacteristicLength(cellgeo);
232
233  // setup initial guess
234  lcoords_guess.resize( gcoords.giveSize() );
235  lcoords_guess.zero();
236
237  // apply Newton-Raphson to solve the problem
238  for ( int nite = 0; nite < 10; nite++ ) {
239  // compute the residual
240  this->local2global(guess, lcoords_guess, cellgeo);
241  res.beDifferenceOf(gcoords, guess);
242
243  // check for convergence
244  error = res.computeNorm();
245  if ( error < convergence_limit ) {
246  break;
247  }
248
249  // compute the corrections
250  this->giveJacobianMatrixAt(jac, lcoords_guess, cellgeo);
251  jac.solveForRhs(res, delta);
252
253  // update guess
255  }
256  if ( error > convergence_limit ) { // Imperfect, could give false negatives.
259  return false;
260  }
261
266
267  bool inside = true;
268  for ( int i = 0; i < 3; i++ ) {
269  if ( answer(i) < ( 0. - POINT_TOL ) ) {
271  inside = false;
272  } else if ( answer(i) > ( 1. + POINT_TOL ) ) {
274  inside = false;
275  }
276  }
277
279  if ( answer(3) < 0. - POINT_TOL ) {
280  return false;
281  }
282  return inside;
283 }
284
285
286 double
288 {
289  return cellgeo.giveVertexCoordinates(1)->distance( cellgeo.giveVertexCoordinates(2) );
290 }
291
292
293 void
294 FEI3dTetQuad :: giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
295 {
296  FloatMatrix dNduvw, coords;
297  this->evaldNdxi(dNduvw, lcoords, cellgeo);
298  coords.resize( 3, dNduvw.giveNumberOfRows() );
299  for ( int i = 1; i <= dNduvw.giveNumberOfRows(); i++ ) {
300  coords.setColumn(* cellgeo.giveVertexCoordinates(i), i);
301  }
302  jacobianMatrix.beProductOf(coords, dNduvw);
303 }
304
305
306 void
307 FEI3dTetQuad :: edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
308 {
309  double xi = lcoords.at(1);
311  answer(0) = 0.5 * ( xi - 1.0 ) * xi;
312  answer(1) = 0.5 * ( xi + 1.0 ) * xi;
313  answer(2) = 1.0 - xi * xi;
314 }
315
316 void
318  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
319 {
320  IntArray edgeNodes;
321  this->computeLocalEdgeMapping(edgeNodes, iedge);
323  OOFEM_ERROR("Not supported");
324 }
325
326 void
328  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
329 {
330  IntArray edgeNodes;
331  FloatArray N;
332  this->computeLocalEdgeMapping(edgeNodes, iedge);
333  this->edgeEvalN(N, iedge, lcoords, cellgeo);
334
336  for ( int i = 0; i < N.giveSize(); ++i ) {
338  }
339 }
340
341
342 double
344 {
345  IntArray edgeNodes;
346  this->computeLocalEdgeMapping(edgeNodes, iedge);
348  OOFEM_ERROR("Not supported");
349  return -1;
350 }
351
352
353 void
355 {
356  edgeNodes.resize(3);
357
358  if ( iedge == 1 ) { // edge between nodes 1 2
359  edgeNodes(0) = 1;
360  edgeNodes(1) = 2;
361  edgeNodes(2) = 5;
362  } else if ( iedge == 2 ) { // edge between nodes 2 3
363  edgeNodes(0) = 2;
364  edgeNodes(1) = 3;
365  edgeNodes(2) = 6;
366  } else if ( iedge == 3 ) { // edge between nodes 3 1
367  edgeNodes(0) = 3;
368  edgeNodes(1) = 1;
369  edgeNodes(2) = 7;
370  } else if ( iedge == 4 ) { // edge between nodes 1 4
371  edgeNodes(0) = 1;
372  edgeNodes(1) = 4;
373  edgeNodes(2) = 8;
374  } else if ( iedge == 5 ) { // edge between nodes 2 4
375  edgeNodes(0) = 2;
376  edgeNodes(1) = 4;
377  edgeNodes(2) = 9;
378  } else if ( iedge == 6 ) { // edge between nodes 3 4
379  edgeNodes(0) = 3;
380  edgeNodes(1) = 4;
381  edgeNodes(2) = 10;
382  } else {
383  OOFEM_ERROR("wrong edge number (%d)", iedge);
384  }
385 }
386
387 double
389 {
391  OOFEM_ERROR("Not supported");
392  return -1;
393 }
394
395 void
396 FEI3dTetQuad :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
397 {
398  double l1 = lcoords.at(1);
399  double l2 = lcoords.at(2);
400  double l3 = 1. - l1 - l2;
401
403
404  answer.at(1) = ( 2. * l1 - 1. ) * l1;
405  answer.at(2) = ( 2. * l2 - 1. ) * l2;
406  answer.at(3) = ( 2. * l3 - 1. ) * l3;
407  answer.at(4) = 4. * l1 * l2;
408  answer.at(5) = 4. * l2 * l3;
409  answer.at(6) = 4. * l3 * l1;
410 }
411
412 void
414  const FloatArray &lcoords, const FEICellGeometry &cellgeo)
415 {
416  IntArray nodes;
417  FloatArray N;
418  this->computeLocalSurfaceMapping(nodes, isurf);
419  this->surfaceEvalN(N, isurf, lcoords, cellgeo);
420
422  for ( int i = 0; i < N.giveSize(); ++i ) {
424  }
425 }
426
427 void
428 FEI3dTetQuad :: surfaceEvaldNdx(FloatMatrix &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
429 {
430  IntArray snodes;
431  this->computeLocalSurfaceMapping(snodes, isurf);
432
433  FloatArray lcoords_tet(4);
434  lcoords_tet.at(snodes.at(1)) = lcoords.at(1);
435  lcoords_tet.at(snodes.at(2)) = lcoords.at(2);
436  lcoords_tet.at(snodes.at(3)) = 1. - lcoords.at(1) - lcoords.at(2);
437
438  FloatMatrix fullB;
439  this->evaldNdx(fullB, lcoords_tet, cellgeo);
441  for ( int i = 1; i <= snodes.giveSize(); ++i ) {
442  for ( int j = 1; j <= 3; ++j ) {
443  answer.at(i, j) = fullB.at(snodes.at(i), j);
444  }
445  }
446 }
447
448 double
449 FEI3dTetQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
450 {
451  IntArray snodes(3);
452  FloatArray a, b;
453  this->computeLocalSurfaceMapping(snodes, isurf);
454
455  double l1, l2, l3;
456  l1 = lcoords.at(1);
457  l2 = lcoords.at(2);
458  l3 = 1.0 - l1 - l2;
459
460  FloatArray dNdxi(6), dNdeta(6);
461
462  dNdxi(0) = 4.0 * l1 - 1.0;
463  dNdxi(1) = 0.0;
464  dNdxi(2) = -1.0 * ( 4.0 * l3 - 1.0 );
465  dNdxi(3) = 4.0 * l2;
466  dNdxi(4) = -4.0 * l2;
467  dNdxi(5) = 4.0 * l3 - 4.0 * l1;
468
469  dNdeta(0) = 0.0;
470  dNdeta(1) = 4.0 * l2 - 1.0;
471  dNdeta(2) = -1.0 * ( 4.0 * l3 - 1.0 );
472  dNdeta(3) = 4.0 * l1;
473  dNdeta(4) = 4.0 * l3 - 4.0 * l2;
474  dNdeta(5) = -4.0 * l1;
475
476  for ( int i = 0; i < 6; ++i ) {
477  a.add( dNdxi(i), * cellgeo.giveVertexCoordinates( snodes(i) ) );
478  b.add( dNdeta(i), * cellgeo.giveVertexCoordinates( snodes(i) ) );
479  }
482 }
483
484 double
486  const FEICellGeometry &cellgeo)
487 {
488  FloatArray normal;
489  return this->surfaceEvalNormal(normal, isurf, lcoords, cellgeo);
490 }
491
492 void
494 {
495  int aNode = 0, bNode = 0, cNode = 0, dNode = 0, eNode = 0, fNode = 0;
496  surfNodes.resize(6);
497
498  if ( isurf == 1 ) {
499  aNode = 1;
500  bNode = 3;
501  cNode = 2;
502  dNode = 7;
503  eNode = 6;
504  fNode = 5;
505  } else if ( isurf == 2 ) {
506  aNode = 1;
507  bNode = 2;
508  cNode = 4;
509  dNode = 5;
510  eNode = 9;
511  fNode = 8;
512  } else if ( isurf == 3 ) {
513  aNode = 2;
514  bNode = 3;
515  cNode = 4;
516  dNode = 6;
517  eNode = 10;
518  fNode = 9;
519  } else if ( isurf == 4 ) {
520  aNode = 1;
521  bNode = 4;
522  cNode = 3;
523  dNode = 8;
524  eNode = 10;
525  fNode = 7;
526  } else {
527  OOFEM_ERROR("wrong surface number (%d)", isurf);
528  }
529
530  surfNodes.at(1) = aNode;
531  surfNodes.at(2) = bNode;
532  surfNodes.at(3) = cNode;
533  surfNodes.at(4) = dNode;
534  surfNodes.at(5) = eNode;
535  surfNodes.at(6) = fNode;
536 }
537
538 double FEI3dTetQuad :: evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
539 {
540  IntArray fNodes;
541  this->computeLocalSurfaceMapping(fNodes, iEdge);
542
543  const FloatArray &c1 = * cellgeo.giveVertexCoordinates( fNodes.at(1) );
544  const FloatArray &c2 = * cellgeo.giveVertexCoordinates( fNodes.at(2) );
545  const FloatArray &c3 = * cellgeo.giveVertexCoordinates( fNodes.at(3) );
546  const FloatArray &c4 = * cellgeo.giveVertexCoordinates( fNodes.at(4) );
547  const FloatArray &c5 = * cellgeo.giveVertexCoordinates( fNodes.at(5) );
548  const FloatArray &c6 = * cellgeo.giveVertexCoordinates( fNodes.at(6) );
549
550  // Expression derived in Mathematica:
551  return (
552  c1(2) * ( c2(1) * ( -2 * c3(0) - 3 * c4(0) + 5 * c5(0) + 5 * c6(0) ) +
553  c3(1) * ( 2 * c2(0) - 5 * c4(0) - 5 * c5(0) + 3 * c6(0) ) +
554  c4(1) * ( 3 * c2(0) + 5 * c3(0) - 4 * c5(0) - 24 * c6(0) ) +
555  c5(1) * ( -5 * c2(0) + 5 * c3(0) + 4 * c4(0) - 4 * c6(0) ) +
556  c6(1) * ( -5 * c2(0) - 3 * c3(0) + 24 * c4(0) + 4 * c5(0) ) ) +
557  c2(2) * ( c1(1) * ( 2 * c3(0) + 3 * c4(0) - 5 * c5(0) - 5 * c6(0) ) +
558  c3(1) * ( -2 * c1(0) + 5 * c4(0) - 3 * c5(0) + 5 * c6(0) ) +
559  c4(1) * ( -3 * c1(0) - 5 * c3(0) + 24 * c5(0) + 4 * c6(0) ) +
560  c5(1) * ( 5 * c1(0) + 3 * c3(0) - 24 * c4(0) - 4 * c6(0) ) +
561  c6(1) * ( 5 * c1(0) - 5 * c3(0) - 4 * c4(0) + 4 * c5(0) ) ) +
562  c3(2) * ( c1(1) * ( -2 * c2(0) + 5 * c4(0) + 5 * c5(0) - 3 * c6(0) ) +
563  c2(1) * ( 2 * c1(0) - 5 * c4(0) + 3 * c5(0) - 5 * c6(0) ) +
564  c4(1) * ( -5 * c1(0) + 5 * c2(0) - 4 * c5(0) + 4 * c6(0) ) +
565  c5(1) * ( -5 * c1(0) - 3 * c2(0) + 4 * c4(0) + 24 * c6(0) ) +
566  c6(1) * ( 3 * c1(0) + 5 * c2(0) - 4 * c4(0) - 24 * c5(0) ) ) +
567  c4(2) * ( c1(1) * ( -3 * c2(0) - 5 * c3(0) + 4 * c5(0) + 24 * c6(0) ) +
568  c2(1) * ( 3 * c1(0) + 5 * c3(0) - 24 * c5(0) - 4 * c6(0) ) +
569  c3(1) * ( 5 * c1(0) - 5 * c2(0) + 4 * c5(0) - 4 * c6(0) ) +
570  c5(1) * ( -4 * c1(0) + 24 * c2(0) - 4 * c3(0) - 16 * c6(0) ) +
571  c6(1) * ( -24 * c1(0) + 4 * c2(0) + 4 * c3(0) + 16 * c5(0) ) ) +
572  c5(2) * ( c1(1) * ( 5 * c2(0) - 5 * c3(0) - 4 * c4(0) + 4 * c6(0) ) +
573  c2(1) * ( -5 * c1(0) - 3 * c3(0) + 24 * c4(0) + 4 * c6(0) ) +
574  c3(1) * ( 5 * c1(0) + 3 * c2(0) - 4 * c4(0) - 24 * c6(0) ) +
575  c4(1) * ( 4 * c1(0) - 24 * c2(0) + 4 * c3(0) + 16 * c6(0) ) +
576  c6(1) * ( -4 * c1(0) - 4 * c2(0) + 24 * c3(0) - 16 * c4(0) ) ) +
577  c6(2) * ( c1(1) * ( 5 * c2(0) + 3 * c3(0) - 24 * c4(0) - 4 * c5(0) ) +
578  c2(1) * ( -5 * c1(0) + 5 * c3(0) + 4 * c4(0) - 4 * c5(0) ) +
579  c3(1) * ( -3 * c1(0) - 5 * c2(0) + 4 * c4(0) + 24 * c5(0) ) +
580  c4(1) * ( 24 * c1(0) - 4 * c2(0) - 4 * c3(0) - 16 * c5(0) ) +
581  c5(1) * ( 4 * c1(0) + 4 * c2(0) - 24 * c3(0) + 16 * c4(0) ) )
582  ) / 30.;
583 }
584
587 {
588  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
589  int points = iRule->getRequiredNumberOfIntegrationPoints(_Tetrahedra, order + 3);
590  iRule->SetUpPointsOnTetrahedra(points, _Unknown);
591  return iRule;
592 }
593
596 {
597  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
598  int points = iRule->getRequiredNumberOfIntegrationPoints(_Triangle, order + 2);
599  iRule->SetUpPointsOnTriangle(points, _Unknown);
600  return iRule;
601 }
602 } // end namespace oofem
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
void beVectorProductOf(const FloatArray &v1, const FloatArray &v2)
Computes vector product (or cross product) of vectors given as parameters, , and stores the result in...
Definition: floatarray.C:415
virtual void evaldNdxi(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
bool solveForRhs(const FloatArray &b, FloatArray &answer, bool transpose=false)
Solves the system of linear equations .
Definition: floatmatrix.C:1112
virtual IntegrationRule * giveBoundaryIntegrationRule(int order, int boundary)
Sets up a suitable integration rule for integrating over the requested boundary.
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual const FloatArray * giveVertexCoordinates(int i) const =0
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
void clear()
Definition: floatarray.h:206
virtual int global2local(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates local coordinates from given global ones.
virtual int SetUpPointsOnTetrahedra(int, MaterialMode mode)
Sets up receiver&#39;s integration points on tetrahedra (volume coords) integration domain.
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
virtual double evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of interpolation functions (shape functions) at given point...
virtual void computeLocalEdgeMapping(IntArray &edgeNodes, int iedge)
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
virtual void surfaceEvaldNdx(FloatMatrix &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of edge interpolation functions (shape functions) at given point...
Abstract base class representing integration rule.
void beDifferenceOf(const FloatArray &a, const FloatArray &b)
Sets receiver to be a - b.
Definition: floatarray.C:341
virtual double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal out of the surface at given point.
double distance(const FloatArray &x) const
Computes the distance between position represented by receiver and position given as parameter...
Definition: floatarray.C:489
virtual void surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of edge interpolation functions (shape functions) at given point.
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
#define POINT_TOL
virtual double giveCharacteristicLength(const FEICellGeometry &cellgeo) const
Returns a characteristic length of the geometry, typically a diagonal or edge length.
virtual double giveVolume(const FEICellGeometry &cellgeo) const
Computes the exact volume.
virtual void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Gives the jacobian matrix at the local coordinates.
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
#define N(p, q)
Definition: mdm.C:367
virtual double surfaceGiveTransformationJacobian(int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge jacobian of transformation between local and global coordinates.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
void resize(int n)
Checks size of receiver towards requested bounds.
Definition: intarray.C:124
Class representing vector of real numbers.
Definition: floatarray.h:82
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
virtual void computeLocalSurfaceMapping(IntArray &edgeNodes, int iedge)
virtual void edgeEvalN(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of edge interpolation functions (shape functions) at given point.
virtual double edgeGiveTransformationJacobian(int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the edge jacobian of transformation between local and global coordinates.
double computeNorm() const
Computes the norm (or length) of the vector.
Definition: floatarray.C:840
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
virtual int getRequiredNumberOfIntegrationPoints(integrationDomain dType, int approxOrder)
Abstract service.
double edgeComputeLength(IntArray &edgeNodes, const FEICellGeometry &cellgeo)
void zero()
Definition: floatarray.C:658
void setColumn(const FloatArray &src, int c)
Sets the values of the matrix in specified column.
Definition: floatmatrix.C:648
virtual void surfaceLocal2global(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.
virtual int SetUpPointsOnTriangle(int, MaterialMode mode)
Sets up receiver&#39;s integration points on triangular (area coords) integration domain.
virtual double evalNXIntegral(int iEdge, const FEICellGeometry &cellgeo)
Computes the integral .
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
int giveSize() const
Definition: intarray.h:203
int giveSize() const
Definition: floatarray.h:218
the oofem namespace is to define a context or scope in which all oofem names are defined.
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
double normalize()
Definition: floatarray.C:828
int giveNumberOfRows() const
Returns number of rows of receiver.
Definition: floatmatrix.h:156
virtual void edgeEvaldNdx(FloatMatrix &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the matrix of derivatives of edge interpolation functions (shape functions) at given point...
virtual void edgeLocal2global(FloatArray &answer, int iedge, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates edge global coordinates from given local ones.