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 #include "itkFEMLoadLandmark.h"
20
21 namespace itk
22 {
23 namespace fem
24 {
25
26 // Overload the CreateAnother() method.
CreateAnother() const27 ::itk::LightObject::Pointer LoadLandmark::CreateAnother() const
28 {
29 ::itk::LightObject::Pointer smartPtr;
30 Pointer copyPtr = Self::New();
31
32 // Copy Load Contents
33 copyPtr->m_Eta = this->m_Eta;
34 copyPtr->m_Point = this->m_Point;
35 copyPtr->m_Target = this->m_Target;
36 copyPtr->m_Source = this->m_Source;
37 copyPtr->m_Force = this->m_Force;
38 copyPtr->m_Solution = this->m_Solution;
39 for(auto i : this->m_Element)
40 {
41 copyPtr->AddNextElement( i );
42 }
43 copyPtr->SetGlobalNumber( this->GetGlobalNumber() );
44
45 smartPtr = static_cast<Pointer>(copyPtr);
46
47 return smartPtr;
48 }
49
50 /**
51 * Find the Element to which the LoadLandmark belongs
52 */
GetAssignedElement(Element::ArrayType1::Pointer elements)53 Element::ConstPointer LoadLandmark::GetAssignedElement(Element::ArrayType1::Pointer elements)
54 {
55 int numElements = elements->Size();
56 for( int n = 0; n < numElements; n++ )
57 {
58 Element::Pointer nel = elements->GetElement(n);
59 if( (nel )->GetLocalFromGlobalCoordinates(m_Source, this->m_Point) )
60 {
61 return dynamic_cast<const Element *>(nel.GetPointer());
62 }
63 }
64
65 return nullptr;
66 }
67
68 /**
69 * Find the Element to which the LoadLandmark belongs
70 */
71
AssignToElement(Element::ArrayType::Pointer elements)72 bool LoadLandmark::AssignToElement(Element::ArrayType::Pointer elements)
73 {
74 bool isFound = false;
75
76 // Compute & store the local coordinates of the undeformed point and
77 // the pointer to the element
78 for( Element::ArrayType::const_iterator n = elements->begin();
79 n != elements->end() && !isFound; ++n )
80 {
81 if( ( *n )->GetLocalFromGlobalCoordinates(m_Source, this->m_Point) )
82 {
83 isFound = true;
84 //std::cout << "Found: " << ( *n ) << std::endl;
85 this->m_Element[0] = *n;
86 }
87 }
88
89 return isFound;
90 }
91
AssignToElement(Element::ArrayType1::Pointer elements)92 bool LoadLandmark::AssignToElement(Element::ArrayType1::Pointer elements)
93 {
94 bool isFound = false;
95
96 // Compute & store the local coordinates of the undeformed point and
97 // the pointer of the element
98
99 int numElements = elements->Size();
100 for( int n = 0;
101 n < numElements && !isFound; n++ )
102 {
103 Element::Pointer nel = elements->GetElement(n);
104 if( (nel )->GetLocalFromGlobalCoordinates(m_Source, this->m_Point) )
105 {
106 isFound = true;
107 //std::cout << "Found: " << nel << std::endl;
108 this->m_Element[0] = nel;
109 }
110 }
111
112 return isFound;
113 }
114
SetEta(double e)115 void LoadLandmark::SetEta(double e)
116 {
117 this->m_Eta = e;
118 }
119
GetEta() const120 double LoadLandmark::GetEta() const
121 {
122 return this->m_Eta;
123 }
124
ApplyLoad(Element::ConstPointer element,Element::VectorType & Fe)125 void LoadLandmark::ApplyLoad(Element::ConstPointer element, Element::VectorType & Fe)
126 {
127 const unsigned int NnDOF = element->GetNumberOfDegreesOfFreedomPerNode();
128 const unsigned int Nnodes = element->GetNumberOfNodes();
129
130 Element::VectorType force(NnDOF, 0.0);
131 Element::VectorType disp(NnDOF, 0.0);
132 Element::VectorType new_source(NnDOF, 0.0);
133 Element::VectorType shapeF;
134
135 Fe.set_size( element->GetNumberOfDegreesOfFreedom() );
136 Fe.fill(0.0);
137
138 // Retrieve the local coordinate at which the force acts
139 Element::VectorType pt = this->GetPoint();
140
141 // Retrieve the stored solution
142 Solution::ConstPointer sol = this->GetSolution();
143
144 // Determine the displacement at point pt
145 constexpr unsigned int TotalSolutionIndex = 1;
146 disp = element->InterpolateSolution(pt, ( *sol ), TotalSolutionIndex);
147
148 // Convert the source to global coordinates
149 new_source = this->GetSource() + disp;
150
151 // Calculate the new force
152 this->SetForce(disp);
153 force = ( this->GetTarget() - new_source ) / this->GetEta();
154
155 // std::cout << " disp " << disp << std::endl;
156 // force /= std::sqrt(fmag);
157 new_source = ( this->GetTarget() - new_source );
158 // std::cout << " force = " << force << " distance " <<
159 // new_source.magnitude() << std::endl;
160
161 Element::Float curdist = new_source.magnitude();
162 if( curdist < 1.0 )
163 {
164 force.fill(0.0);
165 }
166 //std::cout << " LM distance " << curdist << std::endl;
167
168 // "Integrate" at the location of the point load
169 shapeF = element->ShapeFunctions(pt);
170 // Calculate the equivalent nodal loads
171 for( unsigned int n = 0; n < Nnodes; n++ )
172 {
173 for( unsigned int d = 0; d < NnDOF; d++ )
174 {
175 Fe[n * NnDOF + d] += shapeF[n] * force[d];
176 }
177 }
178 }
179
PrintSelf(std::ostream & os,Indent indent) const180 void LoadLandmark::PrintSelf(std::ostream& os, Indent indent) const
181 {
182 Superclass::PrintSelf(os, indent);
183 os << indent << "Eta: " << this->m_Eta << std::endl;
184 os << indent << "Source: " << this->m_Source << std::endl;
185 os << indent << "Target: " << this->m_Target << std::endl;
186 os << indent << "Point: " << this->m_Point << std::endl;
187 os << indent << "Force: " << this->m_Force << std::endl;
188 os << indent << "Solution: " << this->m_Solution << std::endl;
189 }
190
191 } // end namespace fem
192 } // end namespace itk
193