1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 #ifndef itkFEMElement3DStrain_hxx
20 #define itkFEMElement3DStrain_hxx
21
22 #include "itkFEMElement3DStrain.h"
23
24 namespace itk
25 {
26 namespace fem
27 {
28 template <typename TBaseClass>
29 Element3DStrain<TBaseClass>
Element3DStrain()30 ::Element3DStrain() : Superclass(), m_mat(nullptr)
31 {
32 }
33
34 // ////////////////////////////////////////////////////////////////////////
35 /**
36 * Methods related to the physics of the problem.
37 */
38
39 template <typename TBaseClass>
40 void Element3DStrain<TBaseClass>
GetStrainDisplacementMatrix(MatrixType & B,const MatrixType & shapeDgl) const41 ::GetStrainDisplacementMatrix(MatrixType & B, const MatrixType & shapeDgl) const
42 {
43 unsigned int p;
44 unsigned int Nn = 3 * this->GetNumberOfNodes();
45
46 B.set_size(6, Nn);
47
48 // Initialize the B matrix to zero - subsequently, only the nonzero
49 // terms will be filled in
50 B.fill(0.0);
51
52 // Copy the shape function derivatives wrt global coordinates
53 // in right position in B matrix.
54 for( unsigned int i = 0; i < Nn; i++ )
55 {
56 p = i / 3;
57
58 switch( i % 3 )
59 {
60 // case 0: /** Columns 1, 4, 7, ..., 22 */
61 // B[0][i] = shapeDgl[0][p];
62 // B[3][i] = shapeDgl[1][p];
63 // B[5][i] = shapeDgl[2][p];
64 // break;
65
66 // case 1: /** Columns 2, 5, 8, ..., 23 */
67 // B[1][i] = shapeDgl[1][p];
68 // B[3][i] = shapeDgl[0][p];
69 // B[4][i] = shapeDgl[2][p];
70 // break;
71
72 // case 2: /** Columns 3, 6, 9, ..., 24 */
73 // B[2][i] = shapeDgl[2][p];
74 // B[4][i] = shapeDgl[1][p];
75 // B[5][i] = shapeDgl[0][p];
76 // break;
77
78 case 0: /** Columns 1, 4, 7, ..., 22 */
79 B[0][i] = shapeDgl[0][p];
80 B[4][i] = shapeDgl[2][p];
81 B[5][i] = shapeDgl[1][p];
82 break;
83
84 case 1: /** Columns 2, 5, 8, ..., 23 */
85 B[1][i] = shapeDgl[1][p];
86 B[3][i] = shapeDgl[2][p];
87 B[5][i] = shapeDgl[0][p];
88 break;
89
90 case 2: /** Columns 3, 6, 9, ..., 24 */
91 B[2][i] = shapeDgl[2][p];
92 B[3][i] = shapeDgl[1][p];
93 B[4][i] = shapeDgl[0][p];
94 break;
95 }
96 }
97 }
98
99 template <typename TBaseClass>
100 void
101 Element3DStrain<TBaseClass>
GetMaterialMatrix(MatrixType & D) const102 ::GetMaterialMatrix(MatrixType & D) const
103 {
104 D.set_size(6, 6);
105 D.fill(0.0);
106
107 /* Material properties matrix */
108 Float fac = ( m_mat->GetThickness() * m_mat->GetYoungsModulus() )
109 / ( ( 1 + m_mat->GetPoissonsRatio() ) * ( 1 - 2 * m_mat->GetPoissonsRatio() ) );
110 /** Set the elements in the top left quadrant */
111 for( int j = 0; j < 3; j++ )
112 {
113 for( int k = 0; k < 3; k++ )
114 {
115 D[j][k] = m_mat->GetPoissonsRatio();
116 }
117 }
118 /** Set the diagonal elements */
119 for( int k = 0; k < 3; k++ )
120 {
121 D[k][k] = 1 - m_mat->GetPoissonsRatio();
122 }
123 for( int k = 3; k < 6; k++ )
124 {
125 D[k][k] = ( 1 - ( 2 * m_mat->GetPoissonsRatio() ) ) * 0.5;
126 }
127
128 /** Multiply by the factor */
129 D = D * fac;
130 }
131
132 template <typename TBaseClass>
133 void
134 Element3DStrain<TBaseClass>
PrintSelf(std::ostream & os,Indent indent) const135 ::PrintSelf(std::ostream& os, Indent indent) const
136 {
137 Superclass::PrintSelf(os, indent);
138 os << indent << "Materials: " << this->m_mat << std::endl;
139 }
140
141 } // end namespace fem
142 } // end namespace itk
143
144 #endif
145