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: 2009-01-29 21:28:16 $
00007   Version:   $Revision: 1.18 $
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 
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     // We must dynamically cast the given load pointer to the
00060     // correct templated load class, which is given as
00061     // template parameter.
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;/* Need to change if the index changes in CrankNicolsonSolver */
00075     typename Solution::ConstPointer   S=l0->GetSolution(); // has current solution state
00076 
00077     // Order of integration
00078     // FIXME: Allow changing the order of integration by setting a 
00079     //        static member within an element base class.
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         //FIXME temporarily using hexahedron shape f for speed
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       // Adjust the size of a force vector returned from the load object so
00139       // that it is equal to the number of DOFs per node. If the Fg returned
00140       // a vector with less dimensions, we add zero elements. If the Fg
00141       // returned a vector with more dimensions, we remove the extra dimensions.
00142       force.fill(0.0);
00143       
00144       force=l0->Fe(gip,gsol);
00145       // Calculate the equivalent nodal loads
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 }} // end namespace itk::fem
00167 
00168 #endif // #ifndef __itkFEMImageMetricLoadImplementation_h
00169 

Generated at Tue Sep 15 02:55:26 2009 for ITK by doxygen 1.5.8 written by Dimitri van Heesch, © 1997-2000