OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
fei3dhexatriquad.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 - 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
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 "fei3dhexatriquad.h"
36 #include "intarray.h"
37 #include "floatarray.h"
38 #include "floatmatrix.h"
39 #include "gaussintegrationrule.h"
40 
41 namespace oofem {
42 void
43 FEI3dHexaTriQuad :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
44 {
45  double u = lcoords.at(1);
46  double v = lcoords.at(2);
47  double w = lcoords.at(3);
48 
49  double a[] = {
50  0.5 * ( u - 1.0 ) * u, 0.5 * ( u + 1.0 ) * u, 1.0 - u * u
51  };
52  double b[] = {
53  0.5 * ( v - 1.0 ) * v, 0.5 * ( v + 1.0 ) * v, 1.0 - v * v
54  };
55  double c[] = {
56  0.5 * ( w - 1.0 ) * w, 0.5 * ( w + 1.0 ) * w, 1.0 - w * w
57  };
58 
59  answer.resize(27);
60 
61  answer.at(5) = a [ 0 ] * b [ 0 ] * c [ 0 ];
62  answer.at(8) = a [ 1 ] * b [ 0 ] * c [ 0 ];
63  answer.at(16) = a [ 2 ] * b [ 0 ] * c [ 0 ];
64 
65  answer.at(6) = a [ 0 ] * b [ 1 ] * c [ 0 ];
66  answer.at(7) = a [ 1 ] * b [ 1 ] * c [ 0 ];
67  answer.at(14) = a [ 2 ] * b [ 1 ] * c [ 0 ];
68 
69  answer.at(13) = a [ 0 ] * b [ 2 ] * c [ 0 ];
70  answer.at(15) = a [ 1 ] * b [ 2 ] * c [ 0 ];
71  answer.at(22) = a [ 2 ] * b [ 2 ] * c [ 0 ];
72 
73  answer.at(1) = a [ 0 ] * b [ 0 ] * c [ 1 ];
74  answer.at(4) = a [ 1 ] * b [ 0 ] * c [ 1 ];
75  answer.at(12) = a [ 2 ] * b [ 0 ] * c [ 1 ];
76 
77  answer.at(2) = a [ 0 ] * b [ 1 ] * c [ 1 ];
78  answer.at(3) = a [ 1 ] * b [ 1 ] * c [ 1 ];
79  answer.at(10) = a [ 2 ] * b [ 1 ] * c [ 1 ];
80 
81  answer.at(9) = a [ 0 ] * b [ 2 ] * c [ 1 ];
82  answer.at(11) = a [ 1 ] * b [ 2 ] * c [ 1 ];
83  answer.at(21) = a [ 2 ] * b [ 2 ] * c [ 1 ];
84 
85  answer.at(17) = a [ 0 ] * b [ 0 ] * c [ 2 ];
86  answer.at(20) = a [ 1 ] * b [ 0 ] * c [ 2 ];
87  answer.at(26) = a [ 2 ] * b [ 0 ] * c [ 2 ];
88 
89  answer.at(18) = a [ 0 ] * b [ 1 ] * c [ 2 ];
90  answer.at(19) = a [ 1 ] * b [ 1 ] * c [ 2 ];
91  answer.at(24) = a [ 2 ] * b [ 1 ] * c [ 2 ];
92 
93  answer.at(23) = a [ 0 ] * b [ 2 ] * c [ 2 ];
94  answer.at(25) = a [ 1 ] * b [ 2 ] * c [ 2 ];
95  answer.at(27) = a [ 2 ] * b [ 2 ] * c [ 2 ];
96 }
97 
98 
99 void
100 FEI3dHexaTriQuad :: surfaceEvalN(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
101 {
102  double u = lcoords.at(1);
103  double v = lcoords.at(2);
104 
105  double a[] = {
106  0.5 * ( u - 1. ) * u, 0.5 * ( u + 1. ) * u, 1. - u * u
107  };
108  double b[] = {
109  0.5 * ( v - 1. ) * v, 0.5 * ( v + 1. ) * v, 1. - v * v
110  };
111 
112  answer.resize(9);
113 
114  answer.at(1) = a [ 0 ] * b [ 0 ];
115  answer.at(2) = a [ 1 ] * b [ 0 ];
116  answer.at(5) = a [ 2 ] * b [ 0 ];
117 
118  answer.at(4) = a [ 0 ] * b [ 1 ];
119  answer.at(3) = a [ 1 ] * b [ 1 ];
120  answer.at(7) = a [ 2 ] * b [ 1 ];
121 
122  answer.at(8) = a [ 0 ] * b [ 2 ];
123  answer.at(6) = a [ 1 ] * b [ 2 ];
124  answer.at(9) = a [ 2 ] * b [ 2 ];
125 }
126 
127 
128 double
129 FEI3dHexaTriQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
130 {
131  IntArray snodes;
132  FloatArray e1, e2, dNdu(9), dNdv(9);
133 
134  double u = lcoords.at(1);
135  double v = lcoords.at(2);
136 
137  double a[] = {
138  0.5 * ( u - 1. ) * u, 0.5 * ( u + 1. ) * u, 1. - u * u
139  };
140  double b[] = {
141  0.5 * ( v - 1. ) * v, 0.5 * ( v + 1. ) * v, 1. - v * v
142  };
143 
144  double da[] = {
145  u - 0.5, u + 0.5, -2. * u
146  };
147  double db[] = {
148  v - 0.5, v + 0.5, -2. * v
149  };
150 
151  dNdu.at(1) = da [ 0 ] * b [ 0 ];
152  dNdu.at(2) = da [ 1 ] * b [ 0 ];
153  dNdu.at(5) = da [ 2 ] * b [ 0 ];
154 
155  dNdu.at(4) = da [ 0 ] * b [ 1 ];
156  dNdu.at(3) = da [ 1 ] * b [ 1 ];
157  dNdu.at(7) = da [ 2 ] * b [ 1 ];
158 
159  dNdu.at(8) = da [ 0 ] * b [ 2 ];
160  dNdu.at(6) = da [ 1 ] * b [ 2 ];
161  dNdu.at(9) = da [ 2 ] * b [ 2 ];
162 
163 
164  dNdv.at(1) = a [ 0 ] * db [ 0 ];
165  dNdv.at(2) = a [ 1 ] * db [ 0 ];
166  dNdv.at(5) = a [ 2 ] * db [ 0 ];
167 
168  dNdv.at(4) = a [ 0 ] * db [ 1 ];
169  dNdv.at(3) = a [ 1 ] * db [ 1 ];
170  dNdv.at(7) = a [ 2 ] * db [ 1 ];
171 
172  dNdv.at(8) = a [ 0 ] * db [ 2 ];
173  dNdv.at(6) = a [ 1 ] * db [ 2 ];
174  dNdv.at(9) = a [ 2 ] * db [ 2 ];
175 
176  this->computeLocalSurfaceMapping(snodes, isurf);
177  for ( int i = 1; i <= 9; ++i ) {
178  e1.add( dNdu.at(i), * cellgeo.giveVertexCoordinates( snodes.at(i) ) );
179  e2.add( dNdv.at(i), * cellgeo.giveVertexCoordinates( snodes.at(i) ) );
180  }
181 
182  answer.beVectorProductOf(e1, e2);
183  return answer.normalize();
184 }
185 
186 
187 void
188 FEI3dHexaTriQuad :: computeLocalSurfaceMapping(IntArray &nodes, int isurf)
189 {
190  if ( isurf == 1 ) {
191  nodes = { 2, 1, 4, 3, 9, 12, 11, 10, 21};
192  } else if ( isurf == 2 ) {
193  nodes = { 5, 6, 7, 8, 13, 14, 15, 16, 22};
194  } else if ( isurf == 3 ) {
195  nodes = { 1, 2, 6, 5, 9, 18, 13, 17, 23};
196  } else if ( isurf == 4 ) {
197  nodes = { 2, 3, 7, 6, 10, 19, 14, 18, 24};
198  } else if ( isurf == 5 ) {
199  nodes = { 3, 4, 8, 7, 11, 20, 15, 19, 25};
200  } else if ( isurf == 6 ) {
201  nodes = { 4, 1, 5, 8, 12, 17, 16, 20, 26};
202  } else {
203  OOFEM_ERROR("wrong surface number (%d)", isurf);
204  }
205 }
206 
207 
208 void
209 FEI3dHexaTriQuad :: giveLocalDerivative(FloatMatrix &dN, const FloatArray &lcoords)
210 {
211  double u, v, w;
212  u = lcoords.at(1);
213  v = lcoords.at(2);
214  w = lcoords.at(3);
215 
216  // Helpers expressions;
217  double a[] = {
218  0.5 * ( u - 1.0 ) * u, 0.5 * ( u + 1.0 ) * u, 1.0 - u * u
219  };
220  double b[] = {
221  0.5 * ( v - 1.0 ) * v, 0.5 * ( v + 1.0 ) * v, 1.0 - v * v
222  };
223  double c[] = {
224  0.5 * ( w - 1.0 ) * w, 0.5 * ( w + 1.0 ) * w, 1.0 - w * w
225  };
226 
227  double da[] = {
228  -0.5 + u, 0.5 + u, -2.0 * u
229  };
230  double db[] = {
231  -0.5 + v, 0.5 + v, -2.0 * v
232  };
233  double dc[] = {
234  -0.5 + w, 0.5 + w, -2.0 * w
235  };
236 
237  dN.resize(27, 3);
238 
239  dN.at(5, 1) = da [ 0 ] * b [ 0 ] * c [ 0 ];
240  dN.at(8, 1) = da [ 1 ] * b [ 0 ] * c [ 0 ];
241  dN.at(16, 1) = da [ 2 ] * b [ 0 ] * c [ 0 ];
242 
243  dN.at(6, 1) = da [ 0 ] * b [ 1 ] * c [ 0 ];
244  dN.at(7, 1) = da [ 1 ] * b [ 1 ] * c [ 0 ];
245  dN.at(14, 1) = da [ 2 ] * b [ 1 ] * c [ 0 ];
246 
247  dN.at(13, 1) = da [ 0 ] * b [ 2 ] * c [ 0 ];
248  dN.at(15, 1) = da [ 1 ] * b [ 2 ] * c [ 0 ];
249  dN.at(22, 1) = da [ 2 ] * b [ 2 ] * c [ 0 ];
250 
251  dN.at(1, 1) = da [ 0 ] * b [ 0 ] * c [ 1 ];
252  dN.at(4, 1) = da [ 1 ] * b [ 0 ] * c [ 1 ];
253  dN.at(12, 1) = da [ 2 ] * b [ 0 ] * c [ 1 ];
254 
255  dN.at(2, 1) = da [ 0 ] * b [ 1 ] * c [ 1 ];
256  dN.at(3, 1) = da [ 1 ] * b [ 1 ] * c [ 1 ];
257  dN.at(10, 1) = da [ 2 ] * b [ 1 ] * c [ 1 ];
258 
259  dN.at(9, 1) = da [ 0 ] * b [ 2 ] * c [ 1 ];
260  dN.at(11, 1) = da [ 1 ] * b [ 2 ] * c [ 1 ];
261  dN.at(21, 1) = da [ 2 ] * b [ 2 ] * c [ 1 ];
262 
263  dN.at(17, 1) = da [ 0 ] * b [ 0 ] * c [ 2 ];
264  dN.at(20, 1) = da [ 1 ] * b [ 0 ] * c [ 2 ];
265  dN.at(26, 1) = da [ 2 ] * b [ 0 ] * c [ 2 ];
266 
267  dN.at(18, 1) = da [ 0 ] * b [ 1 ] * c [ 2 ];
268  dN.at(19, 1) = da [ 1 ] * b [ 1 ] * c [ 2 ];
269  dN.at(24, 1) = da [ 2 ] * b [ 1 ] * c [ 2 ];
270 
271  dN.at(23, 1) = da [ 0 ] * b [ 2 ] * c [ 2 ];
272  dN.at(25, 1) = da [ 1 ] * b [ 2 ] * c [ 2 ];
273  dN.at(27, 1) = da [ 2 ] * b [ 2 ] * c [ 2 ];
274 
275  //
276 
277  dN.at(5, 2) = a [ 0 ] * db [ 0 ] * c [ 0 ];
278  dN.at(8, 2) = a [ 1 ] * db [ 0 ] * c [ 0 ];
279  dN.at(16, 2) = a [ 2 ] * db [ 0 ] * c [ 0 ];
280 
281  dN.at(6, 2) = a [ 0 ] * db [ 1 ] * c [ 0 ];
282  dN.at(7, 2) = a [ 1 ] * db [ 1 ] * c [ 0 ];
283  dN.at(14, 2) = a [ 2 ] * db [ 1 ] * c [ 0 ];
284 
285  dN.at(13, 2) = a [ 0 ] * db [ 2 ] * c [ 0 ];
286  dN.at(15, 2) = a [ 1 ] * db [ 2 ] * c [ 0 ];
287  dN.at(22, 2) = a [ 2 ] * db [ 2 ] * c [ 0 ];
288 
289  dN.at(1, 2) = a [ 0 ] * db [ 0 ] * c [ 1 ];
290  dN.at(4, 2) = a [ 1 ] * db [ 0 ] * c [ 1 ];
291  dN.at(12, 2) = a [ 2 ] * db [ 0 ] * c [ 1 ];
292 
293  dN.at(2, 2) = a [ 0 ] * db [ 1 ] * c [ 1 ];
294  dN.at(3, 2) = a [ 1 ] * db [ 1 ] * c [ 1 ];
295  dN.at(10, 2) = a [ 2 ] * db [ 1 ] * c [ 1 ];
296 
297  dN.at(9, 2) = a [ 0 ] * db [ 2 ] * c [ 1 ];
298  dN.at(11, 2) = a [ 1 ] * db [ 2 ] * c [ 1 ];
299  dN.at(21, 2) = a [ 2 ] * db [ 2 ] * c [ 1 ];
300 
301  dN.at(17, 2) = a [ 0 ] * db [ 0 ] * c [ 2 ];
302  dN.at(20, 2) = a [ 1 ] * db [ 0 ] * c [ 2 ];
303  dN.at(26, 2) = a [ 2 ] * db [ 0 ] * c [ 2 ];
304 
305  dN.at(18, 2) = a [ 0 ] * db [ 1 ] * c [ 2 ];
306  dN.at(19, 2) = a [ 1 ] * db [ 1 ] * c [ 2 ];
307  dN.at(24, 2) = a [ 2 ] * db [ 1 ] * c [ 2 ];
308 
309  dN.at(23, 2) = a [ 0 ] * db [ 2 ] * c [ 2 ];
310  dN.at(25, 2) = a [ 1 ] * db [ 2 ] * c [ 2 ];
311  dN.at(27, 2) = a [ 2 ] * db [ 2 ] * c [ 2 ];
312 
313  //
314 
315  dN.at(5, 3) = a [ 0 ] * b [ 0 ] * dc [ 0 ];
316  dN.at(8, 3) = a [ 1 ] * b [ 0 ] * dc [ 0 ];
317  dN.at(16, 3) = a [ 2 ] * b [ 0 ] * dc [ 0 ];
318 
319  dN.at(6, 3) = a [ 0 ] * b [ 1 ] * dc [ 0 ];
320  dN.at(7, 3) = a [ 1 ] * b [ 1 ] * dc [ 0 ];
321  dN.at(14, 3) = a [ 2 ] * b [ 1 ] * dc [ 0 ];
322 
323  dN.at(13, 3) = a [ 0 ] * b [ 2 ] * dc [ 0 ];
324  dN.at(15, 3) = a [ 1 ] * b [ 2 ] * dc [ 0 ];
325  dN.at(22, 3) = a [ 2 ] * b [ 2 ] * dc [ 0 ];
326 
327  dN.at(1, 3) = a [ 0 ] * b [ 0 ] * dc [ 1 ];
328  dN.at(4, 3) = a [ 1 ] * b [ 0 ] * dc [ 1 ];
329  dN.at(12, 3) = a [ 2 ] * b [ 0 ] * dc [ 1 ];
330 
331  dN.at(2, 3) = a [ 0 ] * b [ 1 ] * dc [ 1 ];
332  dN.at(3, 3) = a [ 1 ] * b [ 1 ] * dc [ 1 ];
333  dN.at(10, 3) = a [ 2 ] * b [ 1 ] * dc [ 1 ];
334 
335  dN.at(9, 3) = a [ 0 ] * b [ 2 ] * dc [ 1 ];
336  dN.at(11, 3) = a [ 1 ] * b [ 2 ] * dc [ 1 ];
337  dN.at(21, 3) = a [ 2 ] * b [ 2 ] * dc [ 1 ];
338 
339  dN.at(17, 3) = a [ 0 ] * b [ 0 ] * dc [ 2 ];
340  dN.at(20, 3) = a [ 1 ] * b [ 0 ] * dc [ 2 ];
341  dN.at(26, 3) = a [ 2 ] * b [ 0 ] * dc [ 2 ];
342 
343  dN.at(18, 3) = a [ 0 ] * b [ 1 ] * dc [ 2 ];
344  dN.at(19, 3) = a [ 1 ] * b [ 1 ] * dc [ 2 ];
345  dN.at(24, 3) = a [ 2 ] * b [ 1 ] * dc [ 2 ];
346 
347  dN.at(23, 3) = a [ 0 ] * b [ 2 ] * dc [ 2 ];
348  dN.at(25, 3) = a [ 1 ] * b [ 2 ] * dc [ 2 ];
349  dN.at(27, 3) = a [ 2 ] * b [ 2 ] * dc [ 2 ];
350 }
351 
352 
353 double
354 FEI3dHexaTriQuad :: evalNXIntegral(int iSurf, const FEICellGeometry &cellgeo)
355 {
356  IntArray fNodes;
357  this->computeLocalSurfaceMapping(fNodes, iSurf);
358 
359  const FloatArray &c1 = * cellgeo.giveVertexCoordinates( fNodes.at(1) );
360  const FloatArray &c2 = * cellgeo.giveVertexCoordinates( fNodes.at(2) );
361  const FloatArray &c3 = * cellgeo.giveVertexCoordinates( fNodes.at(3) );
362  const FloatArray &c4 = * cellgeo.giveVertexCoordinates( fNodes.at(4) );
363  const FloatArray &c5 = * cellgeo.giveVertexCoordinates( fNodes.at(5) );
364  const FloatArray &c6 = * cellgeo.giveVertexCoordinates( fNodes.at(6) );
365  const FloatArray &c7 = * cellgeo.giveVertexCoordinates( fNodes.at(7) );
366  const FloatArray &c8 = * cellgeo.giveVertexCoordinates( fNodes.at(8) );
367  const FloatArray &c9 = * cellgeo.giveVertexCoordinates( fNodes.at(9) );
368 
369  // Generated with Mathematica (rather unwieldy expression, tried to simplify it as good as possible, but it could probably be better)
370  return (
371  c1(2) * (
372  c2(1) * ( -3 * c3(0) - 3 * c4(0) + 18 * c6(0) - 4 * c7(0) + 18 * c8(0) + 24 * c9(0) ) +
373  c3(1) * ( 3 * c2(0) - 3 * c4(0) + 2 * c5(0) + 2 * c6(0) - 2 * c7(0) - 2 * c8(0) ) +
374  c4(1) * ( 3 * c2(0) + 3 * c3(0) - 18 * c5(0) + 4 * c6(0) - 18 * c7(0) - 24 * c9(0) ) +
375  c5(1) * ( -2 * c3(0) + 18 * c4(0) + 12 * c6(0) + 24 * c7(0) - 108 * c8(0) - 144 * c9(0) ) +
376  c6(1) * ( -18 * c2(0) - 2 * c3(0) - 4 * c4(0) - 12 * c5(0) - 4 * c7(0) + 24 * c8(0) + 16 * c9(0) ) +
377  c7(1) * ( 4 * c2(0) + 2 * c3(0) + 18 * c4(0) - 24 * c5(0) + 4 * c6(0) + 12 * c8(0) - 16 * c9(0) ) +
378  c8(1) * ( -18 * c2(0) + 2 * c3(0) + 108 * c5(0) - 24 * c6(0) - 12 * c7(0) + 144 * c9(0) ) +
379  c9(1) * ( -24 * c2(0) + 24 * c4(0) + 144 * c5(0) - 16 * c6(0) + 16 * c7(0) - 144 * c8(0) )
380  ) +
381  c2(2) * (
382  c1(1) * ( 3 * c3(0) + 3 * c4(0) - 18 * c6(0) + 4 * c7(0) - 18 * c8(0) - 24 * c9(0) ) +
383  c3(1) * ( -3 * c1(0) - 3 * c4(0) + 18 * c5(0) + 18 * c7(0) - 4 * c8(0) + 24 * c9(0) ) +
384  c4(1) * ( -3 * c1(0) + 3 * c3(0) - 2 * c5(0) + 2 * c6(0) + 2 * c7(0) - 2 * c8(0) ) +
385  c5(1) * ( -18 * c3(0) + 2 * c4(0) + 108 * c6(0) - 24 * c7(0) - 12 * c8(0) + 144 * c9(0) ) +
386  c6(1) * ( 18 * c1(0) - 2 * c4(0) - 108 * c5(0) + 12 * c7(0) + 24 * c8(0) - 144 * c9(0) ) +
387  c7(1) * ( -4 * c1(0) - 18 * c3(0) - 2 * c4(0) + 24 * c5(0) - 12 * c6(0) - 4 * c8(0) + 16 * c9(0) ) +
388  c8(1) * ( 18 * c1(0) + 4 * c3(0) + 2 * c4(0) + 12 * c5(0) - 24 * c6(0) + 4 * c7(0) - 16 * c9(0) ) +
389  c9(1) * ( 24 * c1(0) - 24 * c3(0) - 144 * c5(0) + 144 * c6(0) - 16 * c7(0) + 16 * c8(0) )
390  ) +
391  c3(2) * (
392  c1(1) * ( -3 * c2(0) + 3 * c4(0) - 2 * c5(0) - 2 * c6(0) + 2 * c7(0) + 2 * c8(0) ) +
393  c2(1) * ( 3 * c1(0) + 3 * c4(0) - 18 * c5(0) - 18 * c7(0) + 4 * c8(0) - 24 * c9(0) ) +
394  c4(1) * ( -3 * c1(0) - 3 * c2(0) - 4 * c5(0) + 18 * c6(0) + 18 * c8(0) + 24 * c9(0) ) +
395  c5(1) * ( 2 * c1(0) + 18 * c2(0) + 4 * c4(0) + 12 * c6(0) - 24 * c7(0) + 4 * c8(0) - 16 * c9(0) ) +
396  c6(1) * ( 2 * c1(0) - 18 * c4(0) - 12 * c5(0) + 108 * c7(0) - 24 * c8(0) + 144 * c9(0) ) +
397  c7(1) * ( -2 * c1(0) + 18 * c2(0) + 24 * c5(0) - 108 * c6(0) + 12 * c8(0) - 144 * c9(0) ) +
398  c8(1) * ( -2 * c1(0) - 4 * c2(0) - 18 * c4(0) - 4 * c5(0) + 24 * c6(0) - 12 * c7(0) + 16 * c9(0) ) +
399  c9(1) * ( 24 * c2(0) - 24 * c4(0) + 16 * c5(0) - 144 * c6(0) + 144 * c7(0) - 16 * c8(0) )
400  ) +
401  c4(2) * (
402  c1(1) * ( -3 * c2(0) - 3 * c3(0) + 18 * c5(0) - 4 * c6(0) + 18 * c7(0) + 24 * c9(0) ) +
403  c2(1) * ( 3 * c1(0) - 3 * c3(0) + 2 * c5(0) - 2 * c6(0) - 2 * c7(0) + 2 * c8(0) ) +
404  c3(1) * ( 3 * c1(0) + 3 * c2(0) + 4 * c5(0) - 18 * c6(0) - 18 * c8(0) - 24 * c9(0) ) +
405  c5(1) * ( -18 * c1(0) - 2 * c2(0) - 4 * c3(0) - 4 * c6(0) + 24 * c7(0) - 12 * c8(0) + 16 * c9(0) ) +
406  c6(1) * ( 4 * c1(0) + 2 * c2(0) + 18 * c3(0) + 4 * c5(0) + 12 * c7(0) - 24 * c8(0) - 16 * c9(0) ) +
407  c7(1) * ( -18 * c1(0) + 2 * c2(0) - 24 * c5(0) - 12 * c6(0) + 108 * c8(0) + 144 * c9(0) ) +
408  c8(1) * ( -2 * c2(0) + 18 * c3(0) + 12 * c5(0) + 24 * c6(0) - 108 * c7(0) - 144 * c9(0) ) +
409  c9(1) * ( -24 * c1(0) + 24 * c3(0) - 16 * c5(0) + 16 * c6(0) - 144 * c7(0) + 144 * c8(0) )
410  ) +
411  c5(2) * (
412  c1(1) * ( 2 * c3(0) - 18 * c4(0) - 12 * c6(0) - 24 * c7(0) + 108 * c8(0) + 144 * c9(0) ) +
413  c2(1) * ( 18 * c3(0) - 2 * c4(0) - 108 * c6(0) + 24 * c7(0) + 12 * c8(0) - 144 * c9(0) ) +
414  c3(1) * ( -2 * c1(0) - 18 * c2(0) - 4 * c4(0) - 12 * c6(0) + 24 * c7(0) - 4 * c8(0) + 16 * c9(0) ) +
415  c4(1) * ( 18 * c1(0) + 2 * c2(0) + 4 * c3(0) + 4 * c6(0) - 24 * c7(0) + 12 * c8(0) - 16 * c9(0) ) +
416  c6(1) * ( 12 * c1(0) + 108 * c2(0) + 12 * c3(0) - 4 * c4(0) + 32 * c7(0) + 32 * c8(0) - 192 * c9(0) ) +
417  c7(1) * ( 24 * c1(0) - 24 * c2(0) - 24 * c3(0) + 24 * c4(0) - 32 * c6(0) + 32 * c8(0) ) +
418  c8(1) * ( -108 * c1(0) - 12 * c2(0) + 4 * c3(0) - 12 * c4(0) - 32 * c6(0) - 32 * c7(0) + 192 * c9(0) ) +
419  c9(1) * ( -144 * c1(0) + 144 * c2(0) - 16 * c3(0) + 16 * c4(0) + 192 * c6(0) - 192 * c8(0) )
420  ) +
421  c6(2) * (
422  c1(1) * ( 18 * c2(0) + 2 * c3(0) + 4 * c4(0) + 12 * c5(0) + 4 * c7(0) - 24 * c8(0) - 16 * c9(0) ) +
423  c2(1) * ( -18 * c1(0) + 2 * c4(0) + 108 * c5(0) - 12 * c7(0) - 24 * c8(0) + 144 * c9(0) ) +
424  c3(1) * ( -2 * c1(0) + 18 * c4(0) + 12 * c5(0) - 108 * c7(0) + 24 * c8(0) - 144 * c9(0) ) +
425  c4(1) * ( -4 * c1(0) - 2 * c2(0) - 18 * c3(0) - 4 * c5(0) - 12 * c7(0) + 24 * c8(0) + 16 * c9(0) ) +
426  c5(1) * ( -12 * c1(0) - 108 * c2(0) - 12 * c3(0) + 4 * c4(0) - 32 * c7(0) - 32 * c8(0) + 192 * c9(0) ) +
427  c8(1) * ( 24 * c1(0) + 24 * c2(0) - 24 * c3(0) - 24 * c4(0) + 32 * c5(0) - 32 * c7(0) ) +
428  c7(1) * ( -4 * c1(0) + 12 * c2(0) + 108 * c3(0) + 12 * c4(0) + 32 * c5(0) + 32 * c8(0) - 192 * c9(0) ) +
429  c9(1) * ( 16 * c1(0) - 144 * c2(0) + 144 * c3(0) - 16 * c4(0) - 192 * c5(0) + 192 * c7(0) )
430  ) +
431  c7(2) * (
432  c1(1) * ( -4 * c2(0) - 2 * c3(0) - 18 * c4(0) + 24 * c5(0) - 4 * c6(0) - 12 * c8(0) + 16 * c9(0) ) +
433  c2(1) * ( 4 * c1(0) + 18 * c3(0) + 2 * c4(0) - 24 * c5(0) + 12 * c6(0) + 4 * c8(0) - 16 * c9(0) ) +
434  c3(1) * ( 2 * c1(0) - 18 * c2(0) - 24 * c5(0) + 108 * c6(0) - 12 * c8(0) + 144 * c9(0) ) +
435  c4(1) * ( 18 * c1(0) - 2 * c2(0) + 24 * c5(0) + 12 * c6(0) - 108 * c8(0) - 144 * c9(0) ) +
436  c5(1) * ( -24 * c1(0) + 24 * c2(0) + 24 * c3(0) - 24 * c4(0) + 32 * c6(0) - 32 * c8(0) ) +
437  c6(1) * ( 4 * c1(0) - 12 * c2(0) - 108 * c3(0) - 12 * c4(0) - 32 * c5(0) - 32 * c8(0) + 192 * c9(0) ) +
438  c8(1) * ( 12 * c1(0) - 4 * c2(0) + 12 * c3(0) + 108 * c4(0) + 32 * c5(0) + 32 * c6(0) - 192 * c9(0) ) +
439  c9(1) * ( -16 * c1(0) + 16 * c2(0) - 144 * c3(0) + 144 * c4(0) - 192 * c6(0) + 192 * c8(0) )
440  ) +
441  c8(2) * (
442  c1(1) * ( 18 * c2(0) - 2 * c3(0) - 108 * c5(0) + 24 * c6(0) + 12 * c7(0) - 144 * c9(0) ) +
443  c2(1) * ( -18 * c1(0) - 4 * c3(0) - 2 * c4(0) - 12 * c5(0) + 24 * c6(0) - 4 * c7(0) + 16 * c9(0) ) +
444  c3(1) * ( 2 * c1(0) + 4 * c2(0) + 18 * c4(0) + 4 * c5(0) - 24 * c6(0) + 12 * c7(0) - 16 * c9(0) ) +
445  c4(1) * ( 2 * c2(0) - 18 * c3(0) - 12 * c5(0) - 24 * c6(0) + 108 * c7(0) + 144 * c9(0) ) +
446  c5(1) * ( 108 * c1(0) + 12 * c2(0) - 4 * c3(0) + 12 * c4(0) + 32 * c6(0) + 32 * c7(0) - 192 * c9(0) ) +
447  c6(1) * ( -24 * c1(0) - 24 * c2(0) + 24 * c3(0) + 24 * c4(0) - 32 * c5(0) + 32 * c7(0) ) +
448  c7(1) * ( -12 * c1(0) + 4 * c2(0) - 12 * c3(0) - 108 * c4(0) - 32 * c5(0) - 32 * c6(0) + 192 * c9(0) ) +
449  c9(1) * ( 144 * c1(0) - 16 * c2(0) + 16 * c3(0) - 144 * c4(0) + 192 * c5(0) - 192 * c7(0) )
450  ) +
451  c9(2) * (
452  c1(1) * ( 24 * c2(0) - 24 * c4(0) - 144 * c5(0) + 16 * c6(0) - 16 * c7(0) + 144 * c8(0) ) +
453  c2(1) * ( -24 * c1(0) + 24 * c3(0) + 144 * c5(0) - 144 * c6(0) + 16 * c7(0) - 16 * c8(0) ) +
454  c3(1) * ( -24 * c2(0) + 24 * c4(0) - 16 * c5(0) + 144 * c6(0) - 144 * c7(0) + 16 * c8(0) ) +
455  c4(1) * ( 24 * c1(0) - 24 * c3(0) + 16 * c5(0) - 16 * c6(0) + 144 * c7(0) - 144 * c8(0) ) +
456  c5(1) * ( 144 * c1(0) - 144 * c2(0) + 16 * c3(0) - 16 * c4(0) - 192 * c6(0) + 192 * c8(0) ) +
457  c6(1) * ( -16 * c1(0) + 144 * c2(0) - 144 * c3(0) + 16 * c4(0) + 192 * c5(0) - 192 * c7(0) ) +
458  c7(1) * ( 16 * c1(0) - 16 * c2(0) + 144 * c3(0) - 144 * c4(0) + 192 * c6(0) - 192 * c8(0) ) +
459  c8(1) * ( -144 * c1(0) + 16 * c2(0) - 16 * c3(0) + 144 * c4(0) - 192 * c5(0) + 192 * c7(0) ) )
460  ) / 300.0;
461 }
462 
463 IntegrationRule *
464 FEI3dHexaTriQuad :: giveIntegrationRule(int order)
465 {
466  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
468  int points = iRule->getRequiredNumberOfIntegrationPoints(_Cube, order + 15);
469  iRule->SetUpPointsOnCube(points, _Unknown);
470  return iRule;
471 }
472 
473 IntegrationRule *
474 FEI3dHexaTriQuad :: giveBoundaryIntegrationRule(int order, int boundary)
475 {
476  IntegrationRule *iRule = new GaussIntegrationRule(1, NULL);
478  int points = iRule->getRequiredNumberOfIntegrationPoints(_Square, order + 6);
479  iRule->SetUpPointsOnSquare(points, _Unknown);
480  return iRule;
481 }
482 } // end namespace oofem
oofem::FloatArray::beVectorProductOf
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
oofem::FEI3dHexaTriQuad::computeLocalSurfaceMapping
virtual void computeLocalSurfaceMapping(IntArray &nodes, int iSurf)
Definition: fei3dhexatriquad.C:188
oofem::FloatArray::at
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
oofem::FEICellGeometry::giveVertexCoordinates
virtual const FloatArray * giveVertexCoordinates(int i) const =0
oofem::FEICellGeometry
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
oofem::IntegrationRule::SetUpPointsOnSquare
virtual int SetUpPointsOnSquare(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
Definition: integrationrule.h:277
oofem::IntArray
Class implementing an array of integers.
Definition: intarray.h:61
oofem::IntArray::at
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
oofem::IntegrationRule
Abstract base class representing integration rule.
Definition: integrationrule.h:94
oofem::FEI3dHexaTriQuad::giveBoundaryIntegrationRule
virtual IntegrationRule * giveBoundaryIntegrationRule(int order, int boundary)
Sets up a suitable integration rule for integrating over the requested boundary.
Definition: fei3dhexatriquad.C:474
oofem::FEI3dHexaTriQuad::evalNXIntegral
virtual double evalNXIntegral(int iSurf, const FEICellGeometry &cellgeo)
Computes the integral .
Definition: fei3dhexatriquad.C:354
oofem::FEI3dHexaTriQuad::giveLocalDerivative
virtual void giveLocalDerivative(FloatMatrix &dN, const FloatArray &lcoords)
Definition: fei3dhexatriquad.C:209
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition: error.h:61
oofem::IntegrationRule::SetUpPointsOnCube
virtual int SetUpPointsOnCube(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit cube integration domain.
Definition: integrationrule.h:303
oofem::FloatMatrix::at
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
oofem::FloatArray
Class representing vector of real numbers.
Definition: floatarray.h:82
oofem::FEI3dHexaTriQuad::giveIntegrationRule
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
Definition: fei3dhexatriquad.C:464
oofem::FloatMatrix
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
oofem::FEI3dHexaTriQuad::surfaceEvalNormal
virtual double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal out of the surface at given point.
Definition: fei3dhexatriquad.C:129
oofem::FloatMatrix::resize
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
oofem::IntegrationRule::getRequiredNumberOfIntegrationPoints
virtual int getRequiredNumberOfIntegrationPoints(integrationDomain dType, int approxOrder)
Abstract service.
Definition: integrationrule.h:215
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
oofem::FEI3dHexaTriQuad::surfaceEvalN
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.
Definition: fei3dhexatriquad.C:100
oofem::FloatArray::normalize
double normalize()
Normalizes receiver.
Definition: floatarray.C:828
oofem::FEI3dHexaTriQuad::evalN
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
Definition: fei3dhexatriquad.C:43
oofem::FloatArray::add
void add(const FloatArray &src)
Adds array src to receiver.
Definition: floatarray.C:156
oofem::GaussIntegrationRule
Class representing Gaussian-quadrature integration rule.
Definition: gaussintegrationrule.h:59
oofem::FloatArray::resize
void resize(int s)
Resizes receiver towards requested size.
Definition: floatarray.C:631
This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011