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 itkFEMElement2DMembrane_hxx
20 #define itkFEMElement2DMembrane_hxx
21 
22 #include "itkFEMElement2DMembrane.h"
23 
24 namespace itk
25 {
26 namespace fem
27 {
28 template <typename TBaseClass>
29 Element2DMembrane<TBaseClass>
Element2DMembrane()30 ::Element2DMembrane() : 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
41 Element2DMembrane<TBaseClass>
GetStrainDisplacementMatrix(MatrixType & B,const MatrixType & shapeDgl) const42 ::GetStrainDisplacementMatrix(MatrixType & B, const MatrixType & shapeDgl) const
43 {
44   unsigned int p;
45   unsigned int Nn = this->GetNumberOfNodes();
46 
47   B.set_size(4, 2 * Nn); // note minor difference from linear elasticity
48   // Copy the shape function derivatives to the B matrix.
49   for( unsigned int i = 0; i < Nn; i++ )
50     {
51     // Compute B index
52     p = i << 1;
53 
54     // Compute B elements
55     B[0][p]   = shapeDgl[0][i];
56     B[0][p + 1] = 0.0;
57 
58     B[1][p]   = 0.0;
59     B[1][p + 1] = shapeDgl[0][i];
60 
61     B[2][p]   = shapeDgl[1][i];
62     B[2][p + 1] = 0.0;
63 
64     B[3][p]   = 0.0;
65     B[3][p + 1] = shapeDgl[1][i];
66     }
67 }
68 
69 template <typename TBaseClass>
70 void
71 Element2DMembrane<TBaseClass>
GetMassMatrix(MatrixType & Me) const72 ::GetMassMatrix(MatrixType & Me) const
73 {
74   // Call the parent's get matrix function
75   Superclass::GetMassMatrix(Me);
76 
77   // Since parent class doesn't have the material properties,
78   // we need to adjust Me matrix here for the density of the element.
79   Me = Me * m_mat->GetDensityHeatProduct();
80 }
81 
82 template <typename TBaseClass>
83 void
84 Element2DMembrane<TBaseClass>
GetMaterialMatrix(MatrixType & D) const85 ::GetMaterialMatrix(MatrixType & D) const
86 {
87   unsigned int d = 4;
88 
89   D.set_size(d, d);
90 
91   D.fill(0.0);
92 
93   // This is the main difference from the linear elasticity problem.
94   /* Material properties matrix.  Simpler than linear elasticity. */
95   Float disot = m_mat->GetYoungsModulus();
96   for( unsigned int i = 0; i < d; i++ )
97     {
98     D[i][i] = disot;
99     }
100 }
101 
102 template <typename TBaseClass>
103 void
104 Element2DMembrane<TBaseClass>
PrintSelf(std::ostream & os,Indent indent) const105 ::PrintSelf(std::ostream& os, Indent indent) const
106 {
107   Superclass::PrintSelf(os, indent);
108   os << indent << "Materials: " << this->m_mat << std::endl;
109 }
110 
111 } // end namespace fem
112 } // end namespace itk
113 
114 #endif
115