00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __itkFEMImageMetricLoadImplementation_h
00019 #define __itkFEMImageMetricLoadImplementation_h
00020
00021 #include "itkFEMImageMetricLoad.h"
00022
00023 #include "itkFEMElement2DC0LinearLineStress.h"
00024 #include "itkFEMElement2DC1Beam.h"
00025 #include "itkFEMElement2DC0LinearTriangularStress.h"
00026 #include "itkFEMElement2DC0LinearQuadrilateralMembrane.h"
00027 #include "itkFEMElement2DC0LinearQuadrilateralMembrane.h"
00028 #include "itkFEMElement3DC0LinearTetrahedronStrain.h"
00029 #include "itkFEMElement3DC0LinearHexahedronStrain.h"
00030
00031 namespace itk {
00032 namespace fem {
00033
00051 template<class TLoadClass>
00052 class ImageMetricLoadImplementation
00053 {
00054 public:
00055
00056 template<class TElementClassConstPointer>
00057 static void ImplementImageMetricLoad(TElementClassConstPointer element, Element::LoadPointer load, Element::VectorType& Fe )
00058 {
00059
00060
00061
00062 typename TLoadClass::Pointer l0=dynamic_cast<TLoadClass*>(&*load);
00063 if ( !l0 ) throw FEMException(__FILE__, __LINE__, "FEM error");
00064
00065 Implementation(static_cast<Element::ConstPointer>(element),l0,Fe);
00066 }
00067
00068 private:
00069
00070 static const bool m_Registered;
00071
00072 static void Implementation(typename Element::ConstPointer element, typename TLoadClass::Pointer l0, typename Element::VectorType& Fe)
00073 {
00074 const unsigned int TotalSolutionIndex=1;
00075 typename Solution::ConstPointer S=l0->GetSolution();
00076
00077
00078
00079
00080 unsigned int order=l0->GetNumberOfIntegrationPoints();
00081
00082 const unsigned int Nip=element->GetNumberOfIntegrationPoints(order);
00083 const unsigned int Ndofs=element->GetNumberOfDegreesOfFreedomPerNode();
00084 const unsigned int Nnodes=element->GetNumberOfNodes();
00085 unsigned int ImageDimension=Ndofs;
00086
00087 Element::VectorType force(Ndofs,0.0),
00088 ip,gip,gsol,force_tmp,shapef;
00089 Element::Float w,detJ;
00090
00091 Fe.set_size(element->GetNumberOfDegreesOfFreedom());
00092 Fe.fill(0.0);
00093 shapef.set_size(Nnodes);
00094 gsol.set_size(Ndofs);
00095 gip.set_size(Ndofs);
00096
00097 for(unsigned int i=0; i<Nip; i++)
00098 {
00099 element->GetIntegrationPointAndWeight(i,ip,w,order);
00100 if (ImageDimension == 3)
00101 {
00102 #define FASTHEX
00103 #ifdef FASTHEX
00104 float r=ip[0]; float s=ip[1]; float t=ip[2];
00105
00106 shapef[0] = (1 - r) * (1 - s) * (1 - t) * 0.125;
00107 shapef[1] = (1 + r) * (1 - s) * (1 - t) * 0.125;
00108 shapef[2] = (1 + r) * (1 + s) * (1 - t) * 0.125;
00109 shapef[3] = (1 - r) * (1 + s) * (1 - t) * 0.125;
00110 shapef[4] = (1 - r) * (1 - s) * (1 + t) * 0.125;
00111 shapef[5] = (1 + r) * (1 - s) * (1 + t) * 0.125;
00112 shapef[6] = (1 + r) * (1 + s) * (1 + t) * 0.125;
00113 shapef[7] = (1 - r) * (1 + s) * (1 + t) * 0.125;
00114 #else
00115 shapef = element->ShapeFunctions(ip);
00116 #endif
00117 }
00118 else if (ImageDimension==2)
00119 {
00120 shapef = element->ShapeFunctions(ip);
00121 }
00122 float solval,posval;
00123 detJ=element->JacobianDeterminant(ip);
00124
00125 for(unsigned int f=0; f<ImageDimension; f++)
00126 {
00127 solval=0.0;
00128 posval=0.0;
00129 for(unsigned int n=0; n<Nnodes; n++)
00130 {
00131 posval += shapef[n]*((element->GetNodeCoordinates(n))[f]);
00132 solval += shapef[n] * S->GetSolutionValue( element->GetNode(n)->GetDegreeOfFreedom(f) , TotalSolutionIndex);
00133 }
00134 gsol[f]=solval;
00135 gip[f]=posval;
00136 }
00137
00138
00139
00140
00141
00142 force.fill(0.0);
00143
00144 force=l0->Fe(gip,gsol);
00145
00146 for(unsigned int n=0; n<Nnodes; n++)
00147 {
00148 for(unsigned int d=0; d<Ndofs; d++)
00149 {
00150 itk::fem::Element::Float temp=shapef[n]*force[d]*w*detJ;
00151 Fe[n*Ndofs+d] += temp;
00152 }
00153 }
00154
00155 }
00156
00157 }
00158
00159 };
00160
00161
00162 template<class TLoadClass>
00163 const bool ImageMetricLoadImplementation<TLoadClass>::m_Registered = false;
00164
00165
00166 }}
00167
00168 #endif // #ifndef __itkFEMImageMetricLoadImplementation_h
00169