TexGen
Elements.cpp
Go to the documentation of this file.
1/*=============================================================================
2TexGen: Geometric textile modeller.
3Copyright (C) 2006 Martin Sherburn
4
5This program is free software; you can redistribute it and/or
6modify it under the terms of the GNU General Public License
7as published by the Free Software Foundation; either version 2
8of the License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13GNU General Public License for more details.
14
15You should have received a copy of the GNU General Public License
16along with this program; if not, write to the Free Software
17Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18=============================================================================*/
19
20#include "PrecompiledHeaders.h"
21#include "Elements.h"
22
23using namespace TexGen;
24
25CElement::CElement(const CMatrix &P, int iOrder)
26: m_P(P)
27, m_iIntegrationOrder(iOrder)
28{
29}
30
32{
33 m_P = P;
34}
35
37{
38 m_iIntegrationOrder = iOrder;
39}
40
42{
43 vector<INTEGRATION_POINT> IntegrationPoints;
44 GetIntegrationPoints(IntegrationPoints);
45 vector<INTEGRATION_POINT>::iterator itIP;
46 KeMatrix.ZeroMatrix();
47 CMatrix B, D, Kei;
48 GetDMatrix(D);
49 for (itIP = IntegrationPoints.begin(); itIP != IntegrationPoints.end(); ++itIP)
50 {
51 GetBMatrix(B, itIP->Position);
52
53 // Kei = B.T * D * B
54 Kei.EqualsTransposeMultiple(B, D * B);
55
56 KeMatrix += Kei * itIP->dWeight;
57 }
58}
59
60void CElement::GetBMatrix(CMatrix& BMatrix, CMatrix& Position)
61{
62 // B = X * A^-1
63 CMatrix A, X, AInv;
64 GetAMatrix(A);
65 GetXMatrix(X, Position);
66 A.GetInverse(AInv);
67 BMatrix.EqualsMultiple(X, AInv);
68}
69
70void CElementTri::GetIntegrationPoints(vector<INTEGRATION_POINT> &IntegrationPoints)
71{
72 IntegrationPoints.clear();
73 const int iOrder = m_iIntegrationOrder;
75 IP.Position.Initialise(1, 3);
76 switch (iOrder)
77 {
78 case 1:
79 IP.Position(0, 0) = IP.Position(0, 1) = IP.Position(0, 2) = 1.0/3.0;
80 IP.dWeight = 1;
81 IntegrationPoints.push_back(IP);
82 break;
83 case 2:
84 IP.dWeight = 1.0/3.0;
85 IP.Position(0, 1) = IP.Position(0, 2) = 0.5;
86 IP.Position(0, 0) = 0;
87 IntegrationPoints.push_back(IP);
88 IP.Position(0, 0) = IP.Position(0, 2) = 0.5;
89 IP.Position(0, 1) = 0;
90 IntegrationPoints.push_back(IP);
91 IP.Position(0, 0) = IP.Position(0, 1) = 0.5;
92 IP.Position(0, 2) = 0;
93 IntegrationPoints.push_back(IP);
94 break;
95 case 3:
96 IP.dWeight = -27.0/48.0;
97 IP.Position(0, 0) = IP.Position(0, 1) = IP.Position(0, 2) = 1.0/3.0;
98 IntegrationPoints.push_back(IP);
99
100 IP.dWeight = 25.0/48.0;
101 IP.Position(0, 1) = IP.Position(0, 2) = 0.2;
102 IP.Position(0, 0) = 0.6;
103 IntegrationPoints.push_back(IP);
104 IP.Position(0, 0) = IP.Position(0, 2) = 0.2;
105 IP.Position(0, 1) = 0.6;
106 IntegrationPoints.push_back(IP);
107 IP.Position(0, 0) = IP.Position(0, 1) = 0.2;
108 IP.Position(0, 2) = 0.6;
109 IntegrationPoints.push_back(IP);
110 break;
111 case 5:
112 IP.dWeight = 0.255;
113 IP.Position(0, 0) = IP.Position(0, 1) = IP.Position(0, 2) = 1.0/3.0;
114 IntegrationPoints.push_back(IP);
115
116 IP.dWeight = 0.13239415;
117 IP.Position(0, 1) = IP.Position(0, 2) = 0.47014206;
118 IP.Position(0, 0) = 0.05961587;
119 IntegrationPoints.push_back(IP);
120 IP.Position(0, 0) = IP.Position(0, 2) = 0.47014206;
121 IP.Position(0, 1) = 0.05961587;
122 IntegrationPoints.push_back(IP);
123 IP.Position(0, 0) = IP.Position(0, 1) = 0.47014206;
124 IP.Position(0, 2) = 0.05961587;
125 IntegrationPoints.push_back(IP);
126
127 IP.dWeight = 0.12593918;
128 IP.Position(0, 1) = IP.Position(0, 2) = 0.10128651;
129 IP.Position(0, 0) = 0.79742699;
130 IntegrationPoints.push_back(IP);
131 IP.Position(0, 0) = IP.Position(0, 2) = 0.10128651;
132 IP.Position(0, 1) = 0.79742699;
133 IntegrationPoints.push_back(IP);
134 IP.Position(0, 0) = IP.Position(0, 1) = 0.10128651;
135 IP.Position(0, 2) = 0.79742699;
136 IntegrationPoints.push_back(IP);
137 break;
138 }
139 vector<INTEGRATION_POINT>::iterator itIP;
140 for (itIP = IntegrationPoints.begin(); itIP != IntegrationPoints.end(); ++itIP)
141 {
142 itIP->Position = itIP->Position * m_P;
143 }
144}
145
147: CElement(P, iOrder)
148{
149}
150
152{
153 CMatrix A(3, 3);
154 int i, j;
155 for (i=0; i<3; ++i)
156 {
157 for (j=0; j<3; ++j)
158 {
159 if (i!=2)
160 A(i, j) = m_P(j, i);
161 else
162 A(i, j) = 1;
163 }
164 }
165 return 0.5*abs(A.GetDeterminant());
166}
167
171
173: CElementTri(P, iOrder)
174, m_E1(1.0)
175, m_E2(1.0)
176{
177}
178
180{
181 vector<CMatrix> Nodes;
182 CMatrix Node(1, 2);
183 CMatrix Center(1, 2);
184 int i, j;
185 for (i=0; i<3; ++i)
186 {
187 Node(0, 0) = m_P(i, 0);
188 Node(0, 1) = m_P(i, 1);
189 Center += Node;
190 Nodes.push_back(Node);
191 }
192 Center /= 3;
193
194 AMatrix.Initialise(10, 10);
195 vector<CMatrix> Rows;
196 Rows.push_back(GetvMatrix(Nodes[0]));
197 Rows.push_back(GetTheta_xMatrix(Nodes[0]));
198 Rows.push_back(GetTheta_yMatrix(Nodes[0]));
199 Rows.push_back(GetvMatrix(Nodes[1]));
200 Rows.push_back(GetTheta_xMatrix(Nodes[1]));
201 Rows.push_back(GetTheta_yMatrix(Nodes[1]));
202 Rows.push_back(GetvMatrix(Nodes[2]));
203 Rows.push_back(GetTheta_xMatrix(Nodes[2]));
204 Rows.push_back(GetTheta_yMatrix(Nodes[2]));
205 Rows.push_back(GetvMatrix(Center));
206 for (i=0; i<10; ++i)
207 {
208 for (j=0; j<10; ++j)
209 {
210 AMatrix(i, j) = Rows[i](0, j);
211 }
212 }
213}
214
216{
217 XMatrix.Initialise(3, 10);
218 CMatrix R1 = GetEpsilon_xMatrix(Position);
219 CMatrix R2 = GetEpsilon_yMatrix(Position);
220 CMatrix R3 = GetEpsilon_xyMatrix(Position);
221 int j;
222 for (j=0; j<10; ++j)
223 {
224 XMatrix(0, j) = R1(0, j);
225 XMatrix(1, j) = R2(0, j);
226 XMatrix(2, j) = 2*R3(0, j); // Using engineering shear strain here, make sure D matrix also uses this
227 }
228}
229
231{
232/*
233 // This is the isotropic case!
234 const double v = 0.2;
235 const double E = m_E;
236 const double t = 1;
237 DMatrix.Initialise(3, 3);
238 DMatrix(0, 0) = 1;
239 DMatrix(1, 1) = 1;
240 DMatrix(1, 0) = v;
241 DMatrix(0, 1) = v;
242 DMatrix(2, 2) = 0.5*(1-v); // Using engineering shear strain here, make sure X matrix also uses this
243 DMatrix *= (E*t*t*t)/(12*(1-v*v));*/
244 const double v12 = 0.0, v21 = 0.0;
245 const double E1 = m_E1, E2 = m_E2;
246 const double G12 = 0.5*(E1+E2)/(2*(1+0.5*(v12+v21)));
247 const double t = 1;
248 DMatrix.Initialise(3, 3);
249 DMatrix(0, 0) = E1/(1-v12*v21);
250 DMatrix(1, 1) = E2/(1-v12*v21);
251 DMatrix(1, 0) = v12*E2/(1-v12*v21);
252 DMatrix(0, 1) = v12*E2/(1-v12*v21);
253 DMatrix(2, 2) = G12; // Using engineering shear strain here, make sure X matrix also uses this
254// DMatrix *= t*t*t/12;
255
256 // Let's rotate this to the fibre direction
258 if (MatOrient)
259 {
260 Normalise(MatOrient);
261 // Transformation matrix to rotate to fibre direction
262 // From "Mechanics of composite materials" by Robert M. Jones, pg 74
263 // section "2.6. Stress-strain relations for a lamina of arbitrary orientation"
264 CMatrix T(3, 3);
265 T(0, 0) = MatOrient.x*MatOrient.x;
266 T(0, 1) = MatOrient.y*MatOrient.y;
267 T(0, 2) = -2*MatOrient.x*MatOrient.y;
268 T(1, 0) = MatOrient.y*MatOrient.y;
269 T(1, 1) = MatOrient.x*MatOrient.x;
270 T(1, 2) = 2*MatOrient.x*MatOrient.y;
271 T(2, 0) = MatOrient.x*MatOrient.y;
272 T(2, 1) = -MatOrient.x*MatOrient.y;
273 T(2, 2) = MatOrient.x*MatOrient.x-MatOrient.y*MatOrient.y;
274 CMatrix TInv;
275 CMatrix RotatedD;
276 T.GetInverse(TInv);
277
278 // RotatedD = T^-1 * D * T^-T
279 RotatedD.EqualsMultipleTranspose(TInv*DMatrix, TInv);
280
281 DMatrix = RotatedD;
282 }
283}
284
286{
287 double x = Position(0, 0), y = Position(0, 1);
288 CMatrix v(1, 10);
289 v(0, 0) = 1;
290 v(0, 1) = x;
291 v(0, 2) = y;
292 v(0, 3) = x*x;
293 v(0, 4) = x*y;
294 v(0, 5) = y*y;
295 v(0, 6) = x*x*x;
296 v(0, 7) = x*x*y;
297 v(0, 8) = x*y*y;
298 v(0, 9) = y*y*y;
299 return v;
300}
301
303{
304 double x = Position(0, 0), y = Position(0, 1);
305 CMatrix Theta_x(1, 10);
306 Theta_x(0, 0) = 0;
307 Theta_x(0, 1) = 1;
308 Theta_x(0, 2) = 0;
309 Theta_x(0, 3) = 2*x;
310 Theta_x(0, 4) = y;
311 Theta_x(0, 5) = 0;
312 Theta_x(0, 6) = 3*x*x;
313 Theta_x(0, 7) = 2*x*y;
314 Theta_x(0, 8) = y*y;
315 Theta_x(0, 9) = 0;
316 return Theta_x;
317}
318
320{
321 double x = Position(0, 0), y = Position(0, 1);
322 CMatrix Theta_y(1, 10);
323 Theta_y(0, 0) = 0;
324 Theta_y(0, 1) = 0;
325 Theta_y(0, 2) = 1;
326 Theta_y(0, 3) = 0;
327 Theta_y(0, 4) = x;
328 Theta_y(0, 5) = 2*y;
329 Theta_y(0, 6) = 0;
330 Theta_y(0, 7) = x*x;
331 Theta_y(0, 8) = 2*x*y;
332 Theta_y(0, 9) = 3*y*y;
333 return Theta_y;
334}
335
337{
338 double x = Position(0, 0), y = Position(0, 1);
339 CMatrix Epsilon_x(1, 10);
340 Epsilon_x(0, 0) = 0;
341 Epsilon_x(0, 1) = 0;
342 Epsilon_x(0, 2) = 0;
343 Epsilon_x(0, 3) = 2;
344 Epsilon_x(0, 4) = 0;
345 Epsilon_x(0, 5) = 0;
346 Epsilon_x(0, 6) = 6*x;
347 Epsilon_x(0, 7) = 2*y;
348 Epsilon_x(0, 8) = 0;
349 Epsilon_x(0, 9) = 0;
350 return Epsilon_x;
351}
352
354{
355 double x = Position(0, 0), y = Position(0, 1);
356 CMatrix Epsilon_y(1, 10);
357 Epsilon_y(0, 0) = 0;
358 Epsilon_y(0, 1) = 0;
359 Epsilon_y(0, 2) = 0;
360 Epsilon_y(0, 3) = 0;
361 Epsilon_y(0, 4) = 0;
362 Epsilon_y(0, 5) = 2;
363 Epsilon_y(0, 6) = 0;
364 Epsilon_y(0, 7) = 0;
365 Epsilon_y(0, 8) = 2*x;
366 Epsilon_y(0, 9) = 6*y;
367 return Epsilon_y;
368}
369
371{
372 double x = Position(0, 0), y = Position(0, 1);
373 CMatrix Epsilon_xy(1, 10);
374 Epsilon_xy(0, 0) = 0;
375 Epsilon_xy(0, 1) = 0;
376 Epsilon_xy(0, 2) = 0;
377 Epsilon_xy(0, 3) = 0;
378 Epsilon_xy(0, 4) = 1;
379 Epsilon_xy(0, 5) = 0;
380 Epsilon_xy(0, 6) = 0;
381 Epsilon_xy(0, 7) = 2*x;
382 Epsilon_xy(0, 8) = 2*y;
383 Epsilon_xy(0, 9) = 0;
384 return Epsilon_xy;
385}
386
390
392: CElementTri(P, iOrder)
393, m_E(1.0)
394{
395}
396
398{
399 vector<CMatrix> Nodes;
400 CMatrix Node(1, 2);
401 int i, j;
402 for (i=0; i<3; ++i)
403 {
404 Node(0, 0) = m_P(i, 0);
405 Node(0, 1) = m_P(i, 1);
406 Nodes.push_back(Node);
407 }
408
409 AMatrix.Initialise(3, 3);
410 vector<CMatrix> Rows;
411 Rows.push_back(GetvMatrix(Nodes[0]));
412 Rows.push_back(GetvMatrix(Nodes[1]));
413 Rows.push_back(GetvMatrix(Nodes[2]));
414 for (i=0; i<3; ++i)
415 {
416 for (j=0; j<3; ++j)
417 {
418 AMatrix(i, j) = Rows[i](0, j);
419 }
420 }
421}
422
424{
425 XMatrix.Initialise(2, 3);
426 CMatrix R1 = GetEpsilon_xMatrix(Position);
427 CMatrix R2 = GetEpsilon_yMatrix(Position);
428 int j;
429 for (j=0; j<3; ++j)
430 {
431 XMatrix(0, j) = R1(0, j);
432 XMatrix(1, j) = R2(0, j);
433 }
434}
435
437{
438 const double E = m_E;
439 DMatrix.Initialise(2, 2);
440 DMatrix(0, 0) = m_E;
441 DMatrix(1, 1) = 0;
442
443 // Let's rotate this to the fibre direction
445 if (MatOrient)
446 {
447 Normalise(MatOrient);
448 // Transformation matrix to rotate to fibre direction
449 CMatrix T(2, 2);
450 T(0, 0) = MatOrient.x;
451 T(0, 1) = MatOrient.y;
452 T(1, 0) = -MatOrient.y;
453 T(1, 1) = MatOrient.x;
454 CMatrix TInv;
455 CMatrix RotatedD;
456 T.GetInverse(TInv);
457
458 // RotatedD = T^-1 * D * T
459 RotatedD.EqualsMultiple(TInv, DMatrix * T);
460
461 DMatrix = RotatedD;
462 }
463}
464
466{
467 double x = Position(0, 0), y = Position(0, 1);
468 CMatrix v(1, 3);
469 v(0, 0) = 1;
470 v(0, 1) = x;
471 v(0, 2) = y;
472 return v;
473}
474
476{
477 double x = Position(0, 0), y = Position(0, 1);
478 CMatrix v(1, 3);
479 v(0, 0) = 0;
480 v(0, 1) = 1;
481 v(0, 2) = 0;
482 return v;
483}
484
486{
487 double x = Position(0, 0), y = Position(0, 1);
488 CMatrix v(1, 3);
489 v(0, 0) = 0;
490 v(0, 1) = 0;
491 v(0, 2) = 1;
492 return v;
493}
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
Base class for representing finite elements.
Definition: Elements.h:33
int m_iIntegrationOrder
Definition: Elements.h:69
void SetNodeCoordinates(const CMatrix &P)
Definition: Elements.cpp:31
virtual void GetDMatrix(CMatrix &DMatrix)=0
This is the material properties matrix (typically involving E and v)
virtual void GetAMatrix(CMatrix &AMatrix)=0
This is the coordinate matrix given nodal coordinates P.
CMatrix m_P
Definition: Elements.h:68
void GetKeMatrix(CMatrix &KeMatrix)
Definition: Elements.cpp:41
void GetBMatrix(CMatrix &BMatrix, CMatrix &Position)
This is the dimension matrix which is composed of the A and X matrices (B=X*A^-1)
Definition: Elements.cpp:60
void SetIntegrationOrder(int iOrder)
Definition: Elements.cpp:36
virtual void GetXMatrix(CMatrix &XMatrix, CMatrix &Position)=0
This is the matrix which defines strain in terms of the coefficients.
XYZ m_FibreDirection
Definition: Elements.h:70
virtual void GetIntegrationPoints(vector< INTEGRATION_POINT > &IntegrationPoints)=0
Defines a series of integration points.
CMatrix GetEpsilon_xyMatrix(CMatrix &Position)
This is d^2v/dxdy.
Definition: Elements.cpp:370
void GetDMatrix(CMatrix &DMatrix)
This is the material properties matrix (typically involving E and v)
Definition: Elements.cpp:230
CMatrix GetTheta_yMatrix(CMatrix &Position)
This is dv/dy.
Definition: Elements.cpp:319
CMatrix GetvMatrix(CMatrix &Position)
This is the matrix of displacement v in terms of coefficients C1, C2 ... C10.
Definition: Elements.cpp:285
void GetAMatrix(CMatrix &AMatrix)
This is the coordinate matrix given nodal coordinates P.
Definition: Elements.cpp:179
void GetXMatrix(CMatrix &XMatrix, CMatrix &Position)
This is the matrix which defines strain in terms of the coefficients.
Definition: Elements.cpp:215
CElementTriBending(const CMatrix &P=CMatrix(), int iOrder=3)
CElementTriBending ///.
Definition: Elements.cpp:172
CMatrix GetEpsilon_xMatrix(CMatrix &Position)
This is d^2v/dx^2.
Definition: Elements.cpp:336
CMatrix GetTheta_xMatrix(CMatrix &Position)
This is dv/dx.
Definition: Elements.cpp:302
CMatrix GetEpsilon_yMatrix(CMatrix &Position)
This is d^2v/dy^2.
Definition: Elements.cpp:353
Base class for representing triangular elements.
Definition: Elements.h:80
void GetIntegrationPoints(vector< INTEGRATION_POINT > &IntegrationPoints)
Defines a series of integration points.
Definition: Elements.cpp:70
CElementTri(const CMatrix &P, int iOrder)
Definition: Elements.cpp:146
CMatrix GetvMatrix(CMatrix &Position)
This is the matrix of displacement v in terms of coefficients C1, C2 and C3.
Definition: Elements.cpp:465
CMatrix GetEpsilon_xMatrix(CMatrix &Position)
This is dv/dx.
Definition: Elements.cpp:475
void GetAMatrix(CMatrix &AMatrix)
This is the coordinate matrix given nodal coordinates P.
Definition: Elements.cpp:397
void GetXMatrix(CMatrix &XMatrix, CMatrix &Position)
This is the matrix which defines strain in terms of the coefficients.
Definition: Elements.cpp:423
CElementTriTension(const CMatrix &P=CMatrix(), int iOrder=3)
CElementTriTension ///.
Definition: Elements.cpp:391
void GetDMatrix(CMatrix &DMatrix)
This is the material properties matrix (typically involving E and v)
Definition: Elements.cpp:436
CMatrix GetEpsilon_yMatrix(CMatrix &Position)
This is dv/dy.
Definition: Elements.cpp:485
Class to represent a matrix and perform various operations on it.
Definition: Matrix.h:33
void Initialise(int iHeight, int iWidth)
Definition: Matrix.h:144
CMatrix & EqualsMultipleTranspose(const CMatrix &LeftMatrix, const CMatrix &RightMatrix)
This function multiplies the transpose of the left matrix with the right matrix.
Definition: Matrix.h:347
double GetInverse(CMatrix &Inverse) const
Definition: Matrix.h:544
double GetDeterminant() const
Definition: Matrix.h:489
CMatrix & EqualsTransposeMultiple(const CMatrix &LeftMatrix, const CMatrix &RightMatrix)
This function multiplies the transpose of the left matrix with the right matrix.
Definition: Matrix.h:320
void ZeroMatrix()
Definition: Matrix.h:215
CMatrix & EqualsMultiple(const CMatrix &LeftMatrix, const CMatrix &RightMatrix)
This function multiplies the left matrix with the right matrix.
Definition: Matrix.h:293
Namespace containing a series of customised math operations not found in the standard c++ library.
WXYZ & Normalise(WXYZ &Quaternion)
Normalise the quaternion and return it.
Definition: mymath.h:609
Struct for representing an integration point.
Definition: Elements.h:49
Struct for representing points in 2D space.
Definition: mymath.h:103
double x
Definition: mymath.h:104
double y
Definition: mymath.h:104
double x
Definition: mymath.h:57
double y
Definition: mymath.h:57