Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkFEMImageMetricLoadImplementation.h

Go to the documentation of this file.
00001 /*========================================================================= 00002 00003 Program: Insight Segmentation & Registration Toolkit 00004 Module: $RCSfile: itkFEMImageMetricLoadImplementation.h,v $ 00005 Language: C++ 00006 Date: $Date: 2003/09/10 14:29:41 $ 00007 Version: $Revision: 1.16 $ 00008 00009 Copyright (c) Insight Software Consortium. All rights reserved. 00010 See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. 00011 00012 This software is distributed WITHOUT ANY WARRANTY; without even 00013 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00014 PURPOSE. See the above copyright notices for more information. 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 00034 00035 00036 00053 template<class TLoadClass> 00054 class ImageMetricLoadImplementation 00055 { 00056 public: 00057 00058 template<class TElementClassConstPointer> 00059 static void ImplementImageMetricLoad(TElementClassConstPointer element, Element::LoadPointer load, Element::VectorType& Fe ) 00060 { 00061 // We must dynamically cast the given load pointer to the 00062 // correct templated load class, which is given as 00063 // template parameter. 00064 typename TLoadClass::Pointer l0=dynamic_cast<TLoadClass*>(&*load); 00065 if ( !l0 ) throw FEMException(__FILE__, __LINE__, "FEM error"); 00066 00067 Implementation(static_cast<Element::ConstPointer>(element),l0,Fe); 00068 } 00069 00070 private: 00071 00072 static const bool registered; 00073 00074 static void Implementation(typename Element::ConstPointer element, typename TLoadClass::Pointer l0, typename Element::VectorType& Fe) 00075 { 00076 const unsigned int TotalSolutionIndex=1;/* Need to change if the index changes in CrankNicolsonSolver */ 00077 typename Solution::ConstPointer S=l0->GetSolution(); // has current solution state 00078 00079 // Order of integration 00080 // FIXME: Allow changing the order of integration by setting a 00081 // static member within an element base class. 00082 unsigned int order=l0->GetNumberOfIntegrationPoints(); 00083 00084 const unsigned int Nip=element->GetNumberOfIntegrationPoints(order); 00085 const unsigned int Ndofs=element->GetNumberOfDegreesOfFreedomPerNode(); 00086 const unsigned int Nnodes=element->GetNumberOfNodes(); 00087 unsigned int ImageDimension=Ndofs; 00088 00089 Element::VectorType force(Ndofs,0.0), 00090 ip,gip,gsol,force_tmp,shapef; 00091 Element::Float w,detJ; 00092 00093 Fe.resize(element->GetNumberOfDegreesOfFreedom()); 00094 Fe.fill(0.0); 00095 shapef.resize(Nnodes); 00096 gsol.resize(Ndofs); 00097 gip.resize(Ndofs); 00098 00099 for(unsigned int i=0; i<Nip; i++) 00100 { 00101 element->GetIntegrationPointAndWeight(i,ip,w,order); 00102 /* gip=element->GetGlobalFromLocalCoordinates(ip); 00103 gsol=element->InterpolateSolution(ip,*S,TotalSolutionIndex); 00104 //std::cout << gsol << std::endl; 00105 shapeF=element->ShapeFunctions(ip); 00106 */ 00107 00108 if (ImageDimension == 3){ 00109 #define FASTHEX 00110 #ifdef FASTHEX 00111 float r=ip[0]; float s=ip[1]; float t=ip[2]; 00112 //FIXME temporarily using hexahedron shape f for speed 00113 shapef[0] = (1 - r) * (1 - s) * (1 - t) * 0.125; 00114 shapef[1] = (1 + r) * (1 - s) * (1 - t) * 0.125; 00115 shapef[2] = (1 + r) * (1 + s) * (1 - t) * 0.125; 00116 shapef[3] = (1 - r) * (1 + s) * (1 - t) * 0.125; 00117 shapef[4] = (1 - r) * (1 - s) * (1 + t) * 0.125; 00118 shapef[5] = (1 + r) * (1 - s) * (1 + t) * 0.125; 00119 shapef[6] = (1 + r) * (1 + s) * (1 + t) * 0.125; 00120 shapef[7] = (1 - r) * (1 + s) * (1 + t) * 0.125; 00121 #else 00122 shapef = element->ShapeFunctions(ip); 00123 #endif 00124 }else if (ImageDimension==2) shapef = element->ShapeFunctions(ip); 00125 00126 float solval,posval; 00127 detJ=element->JacobianDeterminant(ip); 00128 00129 for(unsigned int f=0; f<ImageDimension; f++) 00130 { 00131 solval=0.0; 00132 posval=0.0; 00133 for(unsigned int n=0; n<Nnodes; n++) 00134 { 00135 posval+=shapef[n]*((element->GetNodeCoordinates(n))[f]); 00136 solval+=shapef[n] * S->GetSolutionValue( element->GetNode(n)->GetDegreeOfFreedom(f) , TotalSolutionIndex); 00137 } 00138 gsol[f]=solval; 00139 gip[f]=posval; 00140 } 00141 00142 // Adjust the size of a force vector returned from the load object so 00143 // that it is equal to the number of DOFs per node. If the Fg returned 00144 // a vector with less dimensions, we add zero elements. If the Fg 00145 // returned a vector with more dimensions, we remove the extra dimensions. 00146 force.fill(0.0); 00147 00148 force=l0->Fe(gip,gsol); 00149 // Calculate the equivalent nodal loads 00150 for(unsigned int n=0; n<Nnodes; n++) 00151 { 00152 for(unsigned int d=0; d<Ndofs; d++) 00153 { 00154 itk::fem::Element::Float temp=shapef[n]*force[d]*w*detJ; 00155 Fe[n*Ndofs+d]+=temp; 00156 } 00157 } 00158 00159 } 00160 00161 } 00162 00163 }; 00164 00165 00166 template<class TLoadClass> 00167 const bool ImageMetricLoadImplementation<TLoadClass>::registered = false; 00168 00169 00170 // When the templated load implementation function is instantiated, 00171 // it will automatically be registered with the VisitorDispatcher so 00172 // that it is called as required. 00173 // Instantiating the implementation function will also instantiate the 00174 // corresponding Load class. 00175 /*template<class TLoadClass> 00176 const bool ImageMetricLoadImplementation<TLoadClass>::registered= 00177 VisitorDispatcher<Element2DC0LinearQuadrilateralMembrane,Element::LoadType, Element2DC0LinearQuadrilateralMembrane::LoadImplementationFunctionPointer> 00178 ::RegisterVisitor((TLoadClass*)0, &ImageMetricLoadImplementation<TLoadClass>::ImplementImageMetricLoad); 00179 */ 00180 00181 00182 00183 }} // end namespace itk::fem 00184 00185 #endif // #ifndef __itkFEMImageMetricLoadImplementation_h

Generated at Sat Mar 31 02:18:28 2007 for ITK by doxygen 1.3.8 written by Dimitri van Heesch, © 1997-2000