Orfeo Toolbox  3.16
itkFEMElement3DC0LinearHexahedron.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkFEMElement3DC0LinearHexahedron.cxx,v $
5  Language: C++
6  Date: $Date: 2009-01-29 20:09:11 $
7  Version: $Revision: 1.8 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 
18 // disable debug warnings in MS compiler
19 #ifdef _MSC_VER
20 #pragma warning(disable: 4786)
21 #endif
22 
23 #include "itkFEMElement3DC0LinearHexahedron.h"
24 #include "vnl/vnl_math.h"
25 
26 namespace itk {
27 namespace fem {
28 
29 void
30 Element3DC0LinearHexahedron
31 ::GetIntegrationPointAndWeight(unsigned int i, VectorType& pt, Float& w, unsigned int order) const
32 {
33  // FIXME: range checking
34 
35  // default integration order=2
36  if (order==0) { order=2; }
37 
38  pt.set_size(3);
39  pt[0] = gaussPoint[order][i%order];
40  pt[1] = gaussPoint[order][(i/order)%order];
41  pt[2] = gaussPoint[order][(i/(order*order))];
42 
43  w=gaussWeight[order][i%order]*gaussWeight[order][(i/order)%order]*gaussWeight[order][(i/(order*order))];
44 
45 }
46 
47 unsigned int
48 Element3DC0LinearHexahedron
49 ::GetNumberOfIntegrationPoints(unsigned int order) const
50 {
51  // FIXME: range checking
52 
53  // default integration order=2
54  if (order==0) { order=2; }
55 
56  return order*order*order;
57 }
58 
59 Element3DC0LinearHexahedron::VectorType
60 Element3DC0LinearHexahedron
61 ::ShapeFunctions( const VectorType& pt ) const
62 {
63  /* Linear hexahedral element has eight shape functions */
64  VectorType shapeF(8);
65 
71  /* given local point x=(r,s,t), where -1 <= r,s,t <= 1 and */
72 
74  shapeF[0] = (1 - pt[0]) * (1 - pt[1]) * (1 - pt[2]) * 0.125;
75 
77  shapeF[1] = (1 + pt[0]) * (1 - pt[1]) * (1 - pt[2]) * 0.125;
78 
80  shapeF[2] = (1 + pt[0]) * (1 + pt[1]) * (1 - pt[2]) * 0.125;
81 
83  shapeF[3] = (1 - pt[0]) * (1 + pt[1]) * (1 - pt[2]) * 0.125;
84 
86  shapeF[4] = (1 - pt[0]) * (1 - pt[1]) * (1 + pt[2]) * 0.125;
87 
89  shapeF[5] = (1 + pt[0]) * (1 - pt[1]) * (1 + pt[2]) * 0.125;
90 
92  shapeF[6] = (1 + pt[0]) * (1 + pt[1]) * (1 + pt[2]) * 0.125;
93 
95  shapeF[7] = (1 - pt[0]) * (1 + pt[1]) * (1 + pt[2]) * 0.125;
96 
97  return shapeF;
98 }
99 
100 void
101 Element3DC0LinearHexahedron
102 ::ShapeFunctionDerivatives( const VectorType& pt, MatrixType& shapeD ) const
103 {
105  shapeD.set_size(3,8);
106 
107  // d(N_1) / d(r)
108  shapeD[0][0] = (-1) * (1 - pt[1]) * (1 - pt[2]) * 0.125;
109 
110  // d(N_1) / d(s)
111  shapeD[1][0] = (-1) * (1 - pt[0]) * (1 - pt[2]) * 0.125;
112 
113  // d(N_1) / d(t)
114  shapeD[2][0] = (-1) * (1 - pt[0]) * (1 - pt[1]) * 0.125;
115 
116  // d(N_2) / d(r)
117  shapeD[0][1] = (+1) * (1 - pt[1]) * (1 - pt[2]) * 0.125;
118 
119  // d(N_2) / d(s)
120  shapeD[1][1] = (-1) * (1 + pt[0]) * (1 - pt[2]) * 0.125;
121 
122  // d(N_2) / d(t)
123  shapeD[2][1] = (-1) * (1 + pt[0]) * (1 - pt[1]) * 0.125;
124 
125  // d(N_3) / d(r)
126  shapeD[0][2] = (+1) * (1 + pt[1]) * (1 - pt[2]) * 0.125;
127 
128  // d(N_3) / d(s)
129  shapeD[1][2] = (+1) * (1 + pt[0]) * (1 - pt[2]) * 0.125;
130 
131  // d(N_3) / d(t)
132  shapeD[2][2] = (-1) * (1 + pt[0]) * (1 + pt[1]) * 0.125;
133 
134  // d(N_4) / d(r)
135  shapeD[0][3] = (-1) * (1 + pt[1]) * (1 - pt[2]) * 0.125;
136 
137  // d(N_4) / d(s)
138  shapeD[1][3] = (+1) * (1 - pt[0]) * (1 - pt[2]) * 0.125;
139 
140  // d(N_4) / d(t)
141  shapeD[2][3] = (-1) * (1 - pt[0]) * (1 + pt[1]) * 0.125;
142 
143  // d(N_5) / d(r)
144  shapeD[0][4] = (-1) * (1 - pt[1]) * (1 + pt[2]) * 0.125;
145 
146  // d(N_5) / d(s)
147  shapeD[1][4] = (-1) * (1 - pt[0]) * (1 + pt[2]) * 0.125;
148 
149  // d(N_5) / d(t)
150  shapeD[2][4] = (+1) * (1 - pt[0]) * (1 - pt[1]) * 0.125;
151 
152  // d(N_6) / d(r)
153  shapeD[0][5] = (+1) * (1 - pt[1]) * (1 + pt[2]) * 0.125;
154 
155  // d(N_6) / d(s)
156  shapeD[1][5] = (-1) * (1 + pt[0]) * (1 + pt[2]) * 0.125;
157 
158  // d(N_6) / d(t)
159  shapeD[2][5] = (+1) * (1 + pt[0]) * (1 - pt[1]) * 0.125;
160 
161  // d(N_7) / d(r)
162  shapeD[0][6] = (+1) * (1 + pt[1]) * (1 + pt[2]) * 0.125;
163 
164  // d(N_7) / d(s)
165  shapeD[1][6] = (+1) * (1 + pt[0]) * (1 + pt[2]) * 0.125;
166 
167  // d(N_7) / d(t)
168  shapeD[2][6] = (+1) * (1 + pt[0]) * (1 + pt[1]) * 0.125;
169 
170  // d(N_8) / d(r)
171  shapeD[0][7] = (-1) * (1 + pt[1]) * (1 + pt[2]) * 0.125;
172 
173  // d(N_8) / d(s)
174  shapeD[1][7] = (+1) * (1 - pt[0]) * (1 + pt[2]) * 0.125;
175 
176  // d(N_8) / d(t)
177  shapeD[2][7] = (+1) * (1 - pt[0]) * (1 + pt[1]) * 0.125;
178 
179 }
180 
181 
182 bool
183 Element3DC0LinearHexahedron
184 ::GetLocalFromGlobalCoordinates( const VectorType& globalPt , VectorType& localPt ) const
185 {
186 
187 // Float x1, x2, x3, x4, y1, y2, y3, y4, xce, yce, xb, yb, xcn, ycn,
188 // A, J1, J2, x0, y0, dx, dy, be, bn, ce, cn;
189 
190  localPt=globalPt;
191  localPt.set_size(3);
192  localPt.fill(0.0);
193 
194  // FIXME!
195 
196  // x1 = this->m_node[0]->GetCoordinates()[0]; y1 = this->m_node[0]->GetCoordinates()[1];
197  // x2 = this->m_node[1]->GetCoordinates()[0]; y2 = this->m_node[1]->GetCoordinates()[1];
198  // x3 = this->m_node[2]->GetCoordinates()[0]; y3 = this->m_node[2]->GetCoordinates()[1];
199  // x4 = this->m_node[3]->GetCoordinates()[0]; y4 = this->m_node[3]->GetCoordinates()[1];
200 
201  // xb = x1 - x2 + x3 - x4;
202  // yb = y1 - y2 + y3 - y4;
203 
204  // xce = x1 + x2 - x3 - x4;
205  // yce = y1 + y2 - y3 - y4;
206 
207  // xcn = x1 - x2 - x3 + x4;
208  // ycn = y1 - y2 - y3 + y4;
209 
210  // A = 0.5 * (((x3 - x1) * (y4 - y2)) - ((x4 - x2) * (y3 - y1)));
211  // J1 = ((x3 - x4) * (y1 - y2)) - ((x1 - x2) * (y3 - y4));
212  // J2 = ((x2 - x3) * (y1 - y4)) - ((x1 - x4) * (y2 - y3));
213 
214  // x0 = 0.25 * (x1 + x2 + x3 + x4);
215  // y0 = 0.25 * (y1 + y2 + y3 + y4);
216 
217  // dx = globalPt[0] - x0;
218  // dy = globalPt[1] - y0;
219 
220  // be = A - (dx * yb) + (dy * xb);
221  // bn = -A - (dx * yb) + (dy * xb);
222  // ce = (dx * yce) - (dy * xce);
223  // cn = (dx * ycn) - (dy * xcn);
224 
225  // localPt[0] = (2 * ce) / (-sqrt((be * be) - (2 * J1 * ce)) - be);
226  // localPt[1] = (2 * cn) / ( vcl_sqrt((bn * bn) + (2 * J2 * cn)) - bn);
227 
228  // FIXME
229  bool IsInside=false;
230 
231  return IsInside;
232 }
233 
237 #ifdef FEM_BUILD_VISUALIZATION
238 void
239 Element3DC0LinearHexahedron
240 ::Draw(CDC* pDC, Solution::ConstPointer sol) const
241 {
242 
243  int x1=m_node[0]->GetCoordinates()[0]*DC_Scale;
244  int y1=m_node[0]->GetCoordinates()[1]*DC_Scale;
245  int z1=m_node[0]->GetCoordinates()[2]*DC_Scale;
246 
247  int x2=m_node[1]->GetCoordinates()[0]*DC_Scale;
248  int y2=m_node[1]->GetCoordinates()[1]*DC_Scale;
249  int z2=m_node[1]->GetCoordinates()[2]*DC_Scale;
250 
251  int x3=m_node[2]->GetCoordinates()[0]*DC_Scale;
252  int y3=m_node[2]->GetCoordinates()[1]*DC_Scale;
253  int z3=m_node[2]->GetCoordinates()[2]*DC_Scale;
254 
255  int x4=m_node[3]->GetCoordinates()[0]*DC_Scale;
256  int y4=m_node[3]->GetCoordinates()[1]*DC_Scale;
257  int z4=m_node[3]->GetCoordinates()[2]*DC_Scale;
258 
259  int x5=m_node[4]->GetCoordinates()[0]*DC_Scale;
260  int y5=m_node[4]->GetCoordinates()[1]*DC_Scale;
261  int z5=m_node[4]->GetCoordinates()[2]*DC_Scale;
262 
263  int x6=m_node[5]->GetCoordinates()[0]*DC_Scale;
264  int y6=m_node[5]->GetCoordinates()[1]*DC_Scale;
265  int z6=m_node[5]->GetCoordinates()[2]*DC_Scale;
266 
267  int x7=m_node[6]->GetCoordinates()[0]*DC_Scale;
268  int y7=m_node[6]->GetCoordinates()[1]*DC_Scale;
269  int z7=m_node[6]->GetCoordinates()[2]*DC_Scale;
270 
271  int x8=m_node[7]->GetCoordinates()[0]*DC_Scale;
272  int y8=m_node[7]->GetCoordinates()[1]*DC_Scale;
273  int z8=m_node[7]->GetCoordinates()[2]*DC_Scale;
274 
275  x1 += sol->GetSolutionValue(this->m_node[0]->GetDegreeOfFreedom(0))*DC_Scale;
276  y1 += sol->GetSolutionValue(this->m_node[0]->GetDegreeOfFreedom(1))*DC_Scale;
277  z1 += sol->GetSolutionValue(this->m_node[0]->GetDegreeOfFreedom(2))*DC_Scale;
278 
279  x2 += sol->GetSolutionValue(this->m_node[1]->GetDegreeOfFreedom(0))*DC_Scale;
280  y2 += sol->GetSolutionValue(this->m_node[1]->GetDegreeOfFreedom(1))*DC_Scale;
281  z2 += sol->GetSolutionValue(this->m_node[1]->GetDegreeOfFreedom(2))*DC_Scale;
282 
283  x3 += sol->GetSolutionValue(this->m_node[2]->GetDegreeOfFreedom(0))*DC_Scale;
284  y3 += sol->GetSolutionValue(this->m_node[2]->GetDegreeOfFreedom(1))*DC_Scale;
285  z3 += sol->GetSolutionValue(this->m_node[2]->GetDegreeOfFreedom(2))*DC_Scale;
286 
287  x4 += sol->GetSolutionValue(this->m_node[3]->GetDegreeOfFreedom(0))*DC_Scale;
288  y4 += sol->GetSolutionValue(this->m_node[3]->GetDegreeOfFreedom(1))*DC_Scale;
289  z4 += sol->GetSolutionValue(this->m_node[3]->GetDegreeOfFreedom(2))*DC_Scale;
290 
291  x5 += sol->GetSolutionValue(this->m_node[4]->GetDegreeOfFreedom(0))*DC_Scale;
292  y5 += sol->GetSolutionValue(this->m_node[4]->GetDegreeOfFreedom(1))*DC_Scale;
293  z5 += sol->GetSolutionValue(this->m_node[4]->GetDegreeOfFreedom(2))*DC_Scale;
294 
295  x6 += sol->GetSolutionValue(this->m_node[5]->GetDegreeOfFreedom(0))*DC_Scale;
296  y6 += sol->GetSolutionValue(this->m_node[5]->GetDegreeOfFreedom(1))*DC_Scale;
297  z6 += sol->GetSolutionValue(this->m_node[5]->GetDegreeOfFreedom(2))*DC_Scale;
298 
299  x7 += sol->GetSolutionValue(this->m_node[6]->GetDegreeOfFreedom(0))*DC_Scale;
300  y7 += sol->GetSolutionValue(this->m_node[6]->GetDegreeOfFreedom(1))*DC_Scale;
301  z7 += sol->GetSolutionValue(this->m_node[6]->GetDegreeOfFreedom(2))*DC_Scale;
302 
303  x8 += sol->GetSolutionValue(this->m_node[7]->GetDegreeOfFreedom(0))*DC_Scale;
304  y8 += sol->GetSolutionValue(this->m_node[7]->GetDegreeOfFreedom(1))*DC_Scale;
305  z8 += sol->GetSolutionValue(this->m_node[7]->GetDegreeOfFreedom(2))*DC_Scale;
306 
307 }
308 #endif
309 
310 }} // end namespace itk::fem
Generated at Sat Feb 2 2013 23:36:41 for Orfeo Toolbox with doxygen 1.8.1.1