Orfeo Toolbox  3.20
Go to the documentation of this file.
1 /*=========================================================================
2
3  Program: Insight Segmentation & Registration Toolkit
5  Language: C++
6  Date: \$Date: 2009-04-05 10:56:50 \$
7  Version: \$Revision: 1.13 \$
8
11
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15
16 =========================================================================*/
17
18 // disable debug warnings in MS compiler
19 #ifdef _MSC_VER
20 #pragma warning(disable: 4786)
21 #endif
22
24 #include "vnl/vnl_math.h"
25
26 namespace itk {
27 namespace fem {
28
29 void
31 ::GetIntegrationPointAndWeight(unsigned int i, VectorType& pt, Float& w, unsigned int order) const
32 {
33  // FIXME: range checking
34
35  // default integration order
36  if (order==0) { order=DefaultIntegrationOrder; }
37
38  pt.set_size(2);
39
40  pt[0]=gaussPoint[order][i%order];
41  pt[1]=gaussPoint[order][i/order];
42
43  w=gaussWeight[order][i%order]*gaussWeight[order][i/order];
44
45 }
46
47 unsigned int
49 ::GetNumberOfIntegrationPoints(unsigned int order) const
50 {
51  // FIXME: range checking
52
53  // default integration order
54  if (order==0) { order=DefaultIntegrationOrder; }
55
56  return order*order;
57 }
58
61 ::ShapeFunctions( const VectorType& pt ) const
62 {
63  /* Linear quadrilateral element has four shape functions */
64  VectorType shapeF(4);
65
71  /* given local point x=(r,s), where -1 <= r,s <= 1 and */
72
73  /* shape function 1: ((1 - r) * (1 - s)) / 4 (node 1) */
74  shapeF[0] = (1 - pt[0]) * (1 - pt[1]) * .25;
75
76  /* shape function 2: ((1 + r) * (1 - s)) / 4 (node 2) */
77  shapeF[1] = (1 + pt[0]) * (1 - pt[1]) * .25;
78
79  /* shape function 3: ((1 + r) * (1 + s)) / 4 (node 3) */
80  shapeF[2] = (1 + pt[0]) * (1 + pt[1]) * .25;
81
82  /* shape function 1: ((1 - r) * (1 + s)) / 4 (node 4) */
83  shapeF[3] = (1 - pt[0]) * (1 + pt[1]) * .25;
84
85  return shapeF;
86 }
87
88 void
91 {
93  shapeD.set_size(2,4);
94
96  shapeD[0][0] = -(1 - pt[1]) * .25;
97
99  shapeD[1][0] = -(1 - pt[0]) * .25;
100
102  shapeD[0][1] = +(1 - pt[1]) * .25;
103
105  shapeD[1][1] = -(1 + pt[0]) * .25;
106
108  shapeD[0][2] = +(1 + pt[1]) * .25;
109
111  shapeD[1][2] = +(1 + pt[0]) * .25;
112
114  shapeD[0][3] = -(1 + pt[1]) * .25;
115
117  shapeD[1][3] = +(1 - pt[0]) * .25;
118
119 }
120
121 bool
123 ::GetLocalFromGlobalCoordinates( const VectorType& globalPt , VectorType& localPt) const
124 {
125
126  Float x1, x2, x3, x4, y1, y2, y3, y4, xce, yce, xb, yb, xcn, ycn,
127  A, J1, J2, x0, y0, dx, dy, be, bn, ce, cn;
128
129  localPt.set_size(2);
130  localPt.fill(0.0);
131
132  x1 = this->m_node[0]->GetCoordinates()[0]; y1 = this->m_node[0]->GetCoordinates()[1];
133  x2 = this->m_node[1]->GetCoordinates()[0]; y2 = this->m_node[1]->GetCoordinates()[1];
134  x3 = this->m_node[2]->GetCoordinates()[0]; y3 = this->m_node[2]->GetCoordinates()[1];
135  x4 = this->m_node[3]->GetCoordinates()[0]; y4 = this->m_node[3]->GetCoordinates()[1];
136
137  xb = x1 - x2 + x3 - x4;
138  yb = y1 - y2 + y3 - y4;
139
140  xce = x1 + x2 - x3 - x4;
141  yce = y1 + y2 - y3 - y4;
142
143  xcn = x1 - x2 - x3 + x4;
144  ycn = y1 - y2 - y3 + y4;
145
146  A = 0.5 * (((x3 - x1) * (y4 - y2)) - ((x4 - x2) * (y3 - y1)));
147  J1 = ((x3 - x4) * (y1 - y2)) - ((x1 - x2) * (y3 - y4));
148  J2 = ((x2 - x3) * (y1 - y4)) - ((x1 - x4) * (y2 - y3));
149
150  x0 = 0.25 * (x1 + x2 + x3 + x4);
151  y0 = 0.25 * (y1 + y2 + y3 + y4);
152
153  dx = globalPt[0] - x0;
154  dy = globalPt[1] - y0;
155
156  be = A - (dx * yb) + (dy * xb);
157  bn = -A - (dx * yb) + (dy * xb);
158  ce = (dx * yce) - (dy * xce);
159  cn = (dx * ycn) - (dy * xcn);
160
161  localPt[0] = (2 * ce) / (-vcl_sqrt((be * be) - (2 * J1 * ce)) - be);
162  localPt[1] = (2 * cn) / ( vcl_sqrt((bn * bn) + (2 * J2 * cn)) - bn);
163
164  bool isInside=true;
165
166  if (localPt[0] < -1.0 || localPt[0] > 1.0 || localPt[1] < -1.0 || localPt[1] > 1.0 )
167  {
168  isInside=false;
169  }
170
171  return isInside;
172 }
173
177 #ifdef FEM_BUILD_VISUALIZATION
178 void
180 ::Draw(CDC* pDC, Solution::ConstPointer sol) const
181 {
182
183  int x1=m_node[0]->GetCoordinates()[0]*DC_Scale;
184  int y1=m_node[0]->GetCoordinates()[1]*DC_Scale;
185
186  int x2=m_node[1]->GetCoordinates()[0]*DC_Scale;
187  int y2=m_node[1]->GetCoordinates()[1]*DC_Scale;
188
189  int x3=m_node[2]->GetCoordinates()[0]*DC_Scale;
190  int y3=m_node[2]->GetCoordinates()[1]*DC_Scale;
191
192  int x4=m_node[3]->GetCoordinates()[0]*DC_Scale;
193  int y4=m_node[3]->GetCoordinates()[1]*DC_Scale;
194
195  x1 += sol->GetSolutionValue(this->m_node[0]->GetDegreeOfFreedom(0))*DC_Scale;
196  y1 += sol->GetSolutionValue(this->m_node[0]->GetDegreeOfFreedom(1))*DC_Scale;
197  x2 += sol->GetSolutionValue(this->m_node[1]->GetDegreeOfFreedom(0))*DC_Scale;
198  y2 += sol->GetSolutionValue(this->m_node[1]->GetDegreeOfFreedom(1))*DC_Scale;
199  x3 += sol->GetSolutionValue(this->m_node[2]->GetDegreeOfFreedom(0))*DC_Scale;
200  y3 += sol->GetSolutionValue(this->m_node[2]->GetDegreeOfFreedom(1))*DC_Scale;
201  x4 += sol->GetSolutionValue(this->m_node[3]->GetDegreeOfFreedom(0))*DC_Scale;
202  y4 += sol->GetSolutionValue(this->m_node[3]->GetDegreeOfFreedom(1))*DC_Scale;
203
204  pDC->MoveTo(x1,y1);
205  pDC->LineTo(x2,y2);
206  pDC->LineTo(x3,y3);
207  pDC->LineTo(x4,y4);
208  pDC->LineTo(x1,y1);
209
210 }
211 #endif
212
213 }} // end namespace itk::fem

Generated at Sat Nov 9 2013 14:28:44 for Orfeo Toolbox with doxygen 1.8.3.1