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 
20 #include "itkFEMElement2DC0LinearLineStress.h"
21 #include "itkFEMSpatialObjectWriter.h"
22 
23 
itkFEMObjectTest(int,char * [])24 int itkFEMObjectTest(int, char *[])
25 {
26   //Need to register default FEM object types,
27   //and setup SpatialReader to recognize FEM types
28   //which is all currently done as a HACK in
29   //the initializaiton of the itk::FEMFactoryBase::GetFactory()
30   itk::FEMFactoryBase::GetFactory()->RegisterDefaultTypes();
31 
32 
33   // Testing the fe mesh validity
34   using FEMObjectType = itk::fem::FEMObject<2>;
35   using FEMObjectTypePointer = FEMObjectType::Pointer;
36 
37   FEMObjectTypePointer femObject = FEMObjectType::New();
38 
39   using NodeType = itk::fem::Element::Node;
40 
41   NodeType::Pointer n1;
42 
43   itk::fem::Element::VectorType pt(2);
44 
45   n1 = NodeType::New();
46   pt[0] = 0.0;
47   pt[1] = 0.0;
48   n1->SetCoordinates(pt);
49   n1->SetGlobalNumber(0);
50   femObject->AddNextNode(n1);
51 
52   n1 = NodeType::New();
53   pt[0] = 1500.0;
54   pt[1] = 0.0;
55   n1->SetCoordinates(pt);
56   n1->SetGlobalNumber(1);
57   femObject->AddNextNode(n1);
58 
59   n1 = NodeType::New();
60   pt[0] = 3000.0;
61   pt[1] = 0.0;
62   n1->SetCoordinates(pt);
63   n1->SetGlobalNumber(2);
64   femObject->AddNextNode(n1);
65 
66   n1 = NodeType::New();
67   pt[0] = 3000.0;
68   pt[1] = 3000.0;
69   n1->SetCoordinates(pt);
70   n1->SetGlobalNumber(3);
71   femObject->AddNextNode(n1);
72 
73   n1 = NodeType::New();
74   pt[0] = 0.0;
75   pt[1] = 4500.0;
76   n1->SetCoordinates(pt);
77   n1->SetGlobalNumber(4);
78   femObject->AddNextNode(n1);
79 
80   itk::fem::MaterialLinearElasticity::Pointer m;
81   m = itk::fem::MaterialLinearElasticity::New();
82   m->SetGlobalNumber(0);               /* Global number of the material */
83   m->SetYoungsModulus(200000000000.0); /* Young modulus */
84   m->SetPoissonsRatio(0.3);
85   m->SetCrossSectionalArea(2000.0); /* Crossection area */
86   m->SetMomentOfInertia(1.0);       /* Momemt of inertia */
87   femObject->AddNextMaterial(m);
88 
89   m = itk::fem::MaterialLinearElasticity::New();
90   m->SetGlobalNumber(1);         /* Global number of the material */
91   m->SetYoungsModulus(200000.0); /* Young modulus */
92   m->SetPoissonsRatio(0.3);
93   m->SetCrossSectionalArea(1200.0); /* Crossection area */
94   m->SetMomentOfInertia(1.0);       /* Momemt of inertia */
95   femObject->AddNextMaterial(m);
96 
97   m = itk::fem::MaterialLinearElasticity::New();
98   m->SetGlobalNumber(2);        /* Global number of the material */
99   m->SetYoungsModulus(70000.0); /* Young modulus */
100   m->SetPoissonsRatio(0.3);
101   m->SetCrossSectionalArea(900.0); /* Crossection area */
102   m->SetMomentOfInertia(1.0);      /* Momemt of inertia */
103   femObject->AddNextMaterial(m);
104 
105   itk::fem::Element2DC0LinearLineStress::Pointer e1;
106 
107   e1 = itk::fem::Element2DC0LinearLineStress::New();
108   e1->SetGlobalNumber(0);
109   e1->SetNode( 0, femObject->GetNode(0) );
110   e1->SetNode( 1, femObject->GetNode(1) );
111   if ( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(1).GetPointer() ))
112     {
113     e1->SetMaterial( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(0).GetPointer() ) );
114     }
115   femObject->AddNextElement( e1);
116 
117   e1 = itk::fem::Element2DC0LinearLineStress::New();
118   e1->SetGlobalNumber(1);
119   e1->SetNode( 0, femObject->GetNode(1) );
120   e1->SetNode( 1, femObject->GetNode(2) );
121   if ( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(0).GetPointer() ))
122     {
123     e1->SetMaterial( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(0).GetPointer() ) );
124     }
125   femObject->AddNextElement( e1);
126 
127   e1 = itk::fem::Element2DC0LinearLineStress::New();
128   e1->SetGlobalNumber(2);
129   e1->SetNode( 0, femObject->GetNode(1) );
130   e1->SetNode( 1, femObject->GetNode(3) );
131   if ( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(2).GetPointer() ))
132     {
133     e1->SetMaterial( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(2).GetPointer() ) );
134     }
135   femObject->AddNextElement( e1);
136 
137   e1 = itk::fem::Element2DC0LinearLineStress::New();
138   e1->SetGlobalNumber(3);
139   e1->SetNode( 0, femObject->GetNode(0) );
140   e1->SetNode( 1, femObject->GetNode(4) );
141   if ( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(1).GetPointer() ))
142     {
143     e1->SetMaterial( dynamic_cast<itk::fem::MaterialLinearElasticity *>( femObject->GetMaterial(1).GetPointer() ) );
144     }
145   femObject->AddNextElement( e1);
146 
147   itk::fem::LoadBC::Pointer l1;
148 
149   l1 = itk::fem::LoadBC::New();
150   l1->SetGlobalNumber(0);
151   l1->SetElement( femObject->GetElement(2) );
152   l1->SetDegreeOfFreedom(2);
153   l1->SetValue( vnl_vector<double>(1, 0.0) );
154   femObject->AddNextLoad( l1 );
155 
156   l1 = itk::fem::LoadBC::New();
157   l1->SetGlobalNumber(1);
158   l1->SetElement( femObject->GetElement(2) );
159   l1->SetDegreeOfFreedom(3);
160   l1->SetValue( vnl_vector<double>(1, 0.0) );
161   femObject->AddNextLoad( l1 );
162 
163   l1 = itk::fem::LoadBC::New();
164   l1->SetGlobalNumber(2);
165   l1->SetElement( femObject->GetElement(3) );
166   l1->SetDegreeOfFreedom(2);
167   l1->SetValue( vnl_vector<double>(1, 0.0) );
168   femObject->AddNextLoad( l1 );
169 
170   l1 = itk::fem::LoadBC::New();
171   l1->SetGlobalNumber(3);
172   l1->SetElement( femObject->GetElement(3) );
173   l1->SetDegreeOfFreedom(3);
174   l1->SetValue( vnl_vector<double>(1, 0.0) );
175   femObject->AddNextLoad( l1 );
176 
177   itk::fem::LoadNode::Pointer l2;
178 
179   l2 = itk::fem::LoadNode::New();
180   l2->SetGlobalNumber(4);
181   l2->SetElement( femObject->GetElement(1) );
182   l2->SetNode(0);
183   vnl_vector<double> F(2);
184   F[0] = 0;
185   F[1] = 30000;
186   l2->SetForce(F);
187   femObject->AddNextLoad( l2 );
188 
189   itk::fem::LoadBCMFC::Pointer bcmfc = itk::fem::LoadBCMFC::New();
190   bcmfc->SetGlobalNumber(5);
191   bcmfc->AddLeftHandSideTerm( itk::fem::LoadBCMFC::MFCTerm(femObject->GetElement(0).GetPointer(), 1, 1) );
192   bcmfc->AddLeftHandSideTerm( itk::fem::LoadBCMFC::MFCTerm(femObject->GetElement(1).GetPointer(), 3, -1) );
193   bcmfc->AddRightHandSideTerm(0.0);
194   femObject->AddNextLoad( bcmfc );
195   femObject->FinalizeMesh();
196 
197   std::cout << "Overall Test : [PASSED]" << std::endl;
198   return EXIT_SUCCESS;
199 }
200