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 "intarray.h"
37 #include "floatarray.h"
38 #include "floatmatrix.h"
39 #include "gaussintegrationrule.h"
40
41 namespace oofem {
42 void
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
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
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
184 }
185
186
187 void
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
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
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
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
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
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 computeLocalSurfaceMapping(IntArray &nodes, int iSurf)
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
virtual int SetUpPointsOnSquare(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit square integration domain.
Class implementing an array of integers.
Definition: intarray.h:61
int & at(int i)
Coefficient access function.
Definition: intarray.h:103
Abstract base class representing integration rule.
virtual IntegrationRule * giveBoundaryIntegrationRule(int order, int boundary)
Sets up a suitable integration rule for integrating over the requested boundary.
virtual double evalNXIntegral(int iSurf, const FEICellGeometry &cellgeo)
Computes the integral .
virtual void giveLocalDerivative(FloatMatrix &dN, const FloatArray &lcoords)
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual int SetUpPointsOnCube(int, MaterialMode mode)
Sets up receiver&#39;s integration points on unit cube integration domain.
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
Class representing vector of real numbers.
Definition: floatarray.h:82
virtual IntegrationRule * giveIntegrationRule(int order)
Sets up a suitable integration rule for numerical integrating over volume.
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
virtual double surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the normal out of the surface at given point.
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.
the oofem namespace is to define a context or scope in which all oofem names are defined.
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.
double normalize()
Definition: floatarray.C:828
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.