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/12/15 14:13:20 $
00007   Version:   $Revision: 1.17 $
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.set_size(element->GetNumberOfDegreesOfFreedom());
00094     Fe.fill(0.0);
00095     shapef.set_size(Nnodes);
00096     gsol.set_size(Ndofs);
00097     gip.set_size(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
00186 

Generated at Mon Mar 12 00:17:29 2007 for ITK by doxygen 1.5.1 written by Dimitri van Heesch, © 1997-2000