1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: itkFEMLoadImplementationGenericLandmarkLoad.cxx
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
// disable debug warnings in MS compiler
#ifdef _MSC_VER
#pragma warning(disable: 4786)
#endif
#include "itkFEMLoadImplementationGenericLandmarkLoad.h"
namespace itk {
namespace fem {
/**
* Handles LandmarkLoad on 2D linear quad stress element
*/
void
LoadImplementationGenericLandmarkLoad
::Implementation(Element::ConstPointer element, LoadLandmark::Pointer load, Element::VectorType& Fe)
{
const unsigned int NnDOF=element->GetNumberOfDegreesOfFreedomPerNode();
const unsigned int Nnodes=element->GetNumberOfNodes();
Element::VectorType force( NnDOF, 0.0 );
Element::VectorType disp( NnDOF, 0.0 );
Element::VectorType new_source (NnDOF, 0.0);
Element::VectorType shapeF;
Fe.set_size(element->GetNumberOfDegreesOfFreedom());
Fe.fill(0.0);
// Retrieve the local coordinate at which the force acts
Element::VectorType pt = load->GetPoint();
// Retrieve the stored solution
Solution::ConstPointer sol = load->GetSolution();
// Determine the displacement at point pt
const unsigned int TotalSolutionIndex=1;
disp = element->InterpolateSolution( pt, (*sol), TotalSolutionIndex );
// Convert the source to global coordinates
new_source = load->GetSource() + disp;
// Calculate the new force
load->m_force = disp;
force = (load->m_target-new_source) / load->eta;
// std::cout << " disp " << disp << std::endl;
//force /= vcl_sqrt(fmag);
new_source = (load->GetTarget() - new_source);
// std::cout << " force = " << force << " distance " << new_source.magnitude() << std::endl;
Element::Float curdist = new_source.magnitude();
if ( curdist < 1.0 )
{
force.fill(0.0);
}
std::cout << " LM distance " << curdist << std::endl;
// "Integrate" at the location of the point load
shapeF = element->ShapeFunctions(pt);
// Calculate the equivalent nodal loads
for(unsigned int n=0; n < Nnodes; n++)
{
for(unsigned int d=0; d < NnDOF; d++)
{
Fe[n*NnDOF+d] += shapeF[n] * force[d];
}
}
}
}} // end namespace itk::fem
|
