<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<html><head><meta name="robots" content="noindex">
<meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>ITK: itkFEMImageMetricLoadImplementation.h Source File</title>
<link href="DoxygenStyle.css" rel="stylesheet" type="text/css">
</head><body bgcolor="#ffffff">
<!-- Section customized for INSIGHT : Tue Jul 17 01:02:45 2001 -->
<center>
<a href="index.html" class="qindex">Main Page</a>
<a href="modules.html" class="qindex">Groups</a>
<a href="namespaces.html" class="qindex">Namespace List</a>
<a href="hierarchy.html" class="qindex">Class Hierarchy</a>
<a href="classes.html" class="qindex">Alphabetical List</a>
<a href="annotated.html" class="qindex">Compound List</a>
<a href="files.html" class="qindex">File
List</a>
<a href="namespacemembers.html" class="qindex">Namespace Members</a>
<a href="functions.html" class="qindex">Compound Members</a>
<a href="globals.html" class="qindex">File Members</a>
<a href="pages.html" class="qindex">Concepts</a></center>
<!-- Generated by Doxygen 1.5.9 -->
<div class="navpath"><a class="el" href="dir_02e23d841c198393a6e961e9a0678fff.html">Code</a> » <a class="el" href="dir_2c33022d19e8077162952ea9ed57656d.html">Numerics</a> » <a class="el" href="dir_316e5ab5c8a9cb56a9dce1c9fc454a23.html">FEM</a>
</div>
<div class="contents">
<h1>itkFEMImageMetricLoadImplementation.h</h1><a href="itkFEMImageMetricLoadImplementation_8h.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">/*=========================================================================</span>
<a name="l00002"></a>00002 <span class="comment"></span>
<a name="l00003"></a>00003 <span class="comment"> Program: Insight Segmentation & Registration Toolkit</span>
<a name="l00004"></a>00004 <span class="comment"> Module: $RCSfile: itkFEMImageMetricLoadImplementation.h,v $</span>
<a name="l00005"></a>00005 <span class="comment"> Language: C++</span>
<a name="l00006"></a>00006 <span class="comment"> Date: $Date: 2009-01-29 21:28:16 $</span>
<a name="l00007"></a>00007 <span class="comment"> Version: $Revision: 1.18 $</span>
<a name="l00008"></a>00008 <span class="comment"></span>
<a name="l00009"></a>00009 <span class="comment"> Copyright (c) Insight Software Consortium. All rights reserved.</span>
<a name="l00010"></a>00010 <span class="comment"> See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.</span>
<a name="l00011"></a>00011 <span class="comment"></span>
<a name="l00012"></a>00012 <span class="comment"> This software is distributed WITHOUT ANY WARRANTY; without even </span>
<a name="l00013"></a>00013 <span class="comment"> the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR </span>
<a name="l00014"></a>00014 <span class="comment"> PURPOSE. See the above copyright notices for more information.</span>
<a name="l00015"></a>00015 <span class="comment"></span>
<a name="l00016"></a>00016 <span class="comment">=========================================================================*/</span>
<a name="l00017"></a>00017
<a name="l00018"></a>00018 <span class="preprocessor">#ifndef __itkFEMImageMetricLoadImplementation_h</span>
<a name="l00019"></a>00019 <span class="preprocessor"></span><span class="preprocessor">#define __itkFEMImageMetricLoadImplementation_h</span>
<a name="l00020"></a>00020 <span class="preprocessor"></span>
<a name="l00021"></a>00021 <span class="preprocessor">#include "<a class="code" href="itkFEMImageMetricLoad_8h.html">itkFEMImageMetricLoad.h</a>"</span>
<a name="l00022"></a>00022
<a name="l00023"></a>00023 <span class="preprocessor">#include "<a class="code" href="itkFEMElement2DC0LinearLineStress_8h.html">itkFEMElement2DC0LinearLineStress.h</a>"</span>
<a name="l00024"></a>00024 <span class="preprocessor">#include "<a class="code" href="itkFEMElement2DC1Beam_8h.html">itkFEMElement2DC1Beam.h</a>"</span>
<a name="l00025"></a>00025 <span class="preprocessor">#include "<a class="code" href="itkFEMElement2DC0LinearTriangularStress_8h.html">itkFEMElement2DC0LinearTriangularStress.h</a>"</span>
<a name="l00026"></a>00026 <span class="preprocessor">#include "<a class="code" href="itkFEMElement2DC0LinearQuadrilateralMembrane_8h.html">itkFEMElement2DC0LinearQuadrilateralMembrane.h</a>"</span>
<a name="l00027"></a>00027 <span class="preprocessor">#include "<a class="code" href="itkFEMElement2DC0LinearQuadrilateralMembrane_8h.html">itkFEMElement2DC0LinearQuadrilateralMembrane.h</a>"</span>
<a name="l00028"></a>00028 <span class="preprocessor">#include "<a class="code" href="itkFEMElement3DC0LinearTetrahedronStrain_8h.html">itkFEMElement3DC0LinearTetrahedronStrain.h</a>"</span>
<a name="l00029"></a>00029 <span class="preprocessor">#include "<a class="code" href="itkFEMElement3DC0LinearHexahedronStrain_8h.html">itkFEMElement3DC0LinearHexahedronStrain.h</a>"</span>
<a name="l00030"></a>00030
<a name="l00031"></a>00031 <span class="keyword">namespace </span>itk {
<a name="l00032"></a>00032 <span class="keyword">namespace </span>fem {
<a name="l00033"></a>00033
<a name="l00051"></a>00051 <span class="keyword">template</span><<span class="keyword">class</span> TLoadClass>
<a name="l00052"></a><a class="code" href="classitk_1_1fem_1_1ImageMetricLoadImplementation.html">00052</a> <span class="keyword">class </span><a class="code" href="classitk_1_1fem_1_1ImageMetricLoadImplementation.html">ImageMetricLoadImplementation</a>
<a name="l00053"></a>00053 {
<a name="l00054"></a>00054 <span class="keyword">public</span>:
<a name="l00055"></a>00055
<a name="l00056"></a>00056 <span class="keyword">template</span><<span class="keyword">class</span> TElementClassConstPo<span class="keywordtype">int</span>er>
<a name="l00057"></a><a class="code" href="classitk_1_1fem_1_1ImageMetricLoadImplementation.html#164c833f8f90d35dc009926452321c23">00057</a> <span class="keyword">static</span> <span class="keywordtype">void</span> <a class="code" href="classitk_1_1fem_1_1ImageMetricLoadImplementation.html#164c833f8f90d35dc009926452321c23">ImplementImageMetricLoad</a>(TElementClassConstPointer element, <a class="code" href="classitk_1_1fem_1_1FEMLightObject.html" title="Base class for all classes that define the FEM system.">Element::LoadPointer</a> load, <a class="code" href="classvnl__vector.html">Element::VectorType</a>& Fe )
<a name="l00058"></a>00058 {
<a name="l00059"></a>00059 <span class="comment">// We must dynamically cast the given load pointer to the</span>
<a name="l00060"></a>00060 <span class="comment">// correct templated load class, which is given as</span>
<a name="l00061"></a>00061 <span class="comment">// template parameter.</span>
<a name="l00062"></a>00062 <span class="keyword">typename</span> <a class="code" href="itkFEMMacro_8h.html#539cce1a3282ba59952dedcbf9cdb23f">TLoadClass::Pointer</a> l0=<span class="keyword">dynamic_cast<</span>TLoadClass*<span class="keyword">></span>(&*load);
<a name="l00063"></a>00063 <span class="keywordflow">if</span> ( !l0 ) <span class="keywordflow">throw</span> <a class="code" href="classitk_1_1fem_1_1FEMException.html" title="Base class for all exception&#39;s that can occur within FEM classes.">FEMException</a>(__FILE__, __LINE__, <span class="stringliteral">"FEM error"</span>);
<a name="l00064"></a>00064
<a name="l00065"></a>00065 Implementation(static_cast<Element::ConstPointer>(element),l0,Fe);
<a name="l00066"></a>00066 }
<a name="l00067"></a>00067
<a name="l00068"></a>00068 <span class="keyword">private</span>:
<a name="l00069"></a>00069
<a name="l00070"></a>00070 <span class="keyword">static</span> <span class="keyword">const</span> <span class="keywordtype">bool</span> m_Registered;
<a name="l00071"></a>00071
<a name="l00072"></a>00072 <span class="keyword">static</span> <span class="keywordtype">void</span> Implementation(<span class="keyword">typename</span> <a class="code" href="classitk_1_1fem_1_1Element.html" title="Abstract base element class.">Element::ConstPointer</a> element, <span class="keyword">typename</span> <a class="code" href="itkFEMMacro_8h.html#539cce1a3282ba59952dedcbf9cdb23f">TLoadClass::Pointer</a> l0, <span class="keyword">typename</span> <a class="code" href="classvnl__vector.html">Element::VectorType</a>& Fe)
<a name="l00073"></a>00073 {
<a name="l00074"></a>00074 <span class="keyword">const</span> <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> TotalSolutionIndex=1;<span class="comment">/* Need to change if the index changes in CrankNicolsonSolver */</span>
<a name="l00075"></a>00075 <span class="keyword">typename</span> <a class="code" href="classitk_1_1fem_1_1Solution.html" title="Provides functions to access the values of the solution vector.">Solution::ConstPointer</a> S=l0->GetSolution(); <span class="comment">// has current solution state</span>
<a name="l00076"></a>00076
<a name="l00077"></a>00077 <span class="comment">// Order of integration</span>
<a name="l00078"></a>00078 <span class="comment">// FIXME: Allow changing the order of integration by setting a </span>
<a name="l00079"></a>00079 <span class="comment">// static member within an element base class.</span>
<a name="l00080"></a>00080 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> order=l0->GetNumberOfIntegrationPoints();
<a name="l00081"></a>00081
<a name="l00082"></a>00082 <span class="keyword">const</span> <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> Nip=element-><a class="code" href="classitk_1_1fem_1_1Element.html#84f91e962f1cddd07778d98058658cac">GetNumberOfIntegrationPoints</a>(order);
<a name="l00083"></a>00083 <span class="keyword">const</span> <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> Ndofs=element-><a class="code" href="classitk_1_1fem_1_1Element.html#b45acf8dbc7b93a0582451b6f4decc20">GetNumberOfDegreesOfFreedomPerNode</a>();
<a name="l00084"></a>00084 <span class="keyword">const</span> <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> Nnodes=element-><a class="code" href="classitk_1_1fem_1_1Element.html#a149d5efdebe193798e310898cd0196f">GetNumberOfNodes</a>();
<a name="l00085"></a>00085 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> ImageDimension=Ndofs;
<a name="l00086"></a>00086
<a name="l00087"></a>00087 <a class="code" href="classvnl__vector.html">Element::VectorType</a> force(Ndofs,0.0),
<a name="l00088"></a>00088 ip,gip,gsol,force_tmp,shapef;
<a name="l00089"></a>00089 <a class="code" href="classitk_1_1fem_1_1Element.html#c55f879c15daa302723012b7f9d21b94">Element::Float</a> w,detJ;
<a name="l00090"></a>00090
<a name="l00091"></a>00091 Fe.set_size(element-><a class="code" href="classitk_1_1fem_1_1Element.html#b3199f871a48d5df67247c40182da482">GetNumberOfDegreesOfFreedom</a>());
<a name="l00092"></a>00092 Fe.fill(0.0);
<a name="l00093"></a>00093 shapef.set_size(Nnodes);
<a name="l00094"></a>00094 gsol.set_size(Ndofs);
<a name="l00095"></a>00095 gip.set_size(Ndofs);
<a name="l00096"></a>00096
<a name="l00097"></a>00097 <span class="keywordflow">for</span>(<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i=0; i<Nip; i++)
<a name="l00098"></a>00098 {
<a name="l00099"></a>00099 element-><a class="code" href="classitk_1_1fem_1_1Element.html#237a942ff6fc96724817b3ba0d401bce">GetIntegrationPointAndWeight</a>(i,ip,w,order);
<a name="l00100"></a>00100 <span class="keywordflow">if</span> (ImageDimension == 3)
<a name="l00101"></a>00101 {
<a name="l00102"></a>00102 <span class="preprocessor">#define FASTHEX</span>
<a name="l00103"></a>00103 <span class="preprocessor"></span><span class="preprocessor">#ifdef FASTHEX</span>
<a name="l00104"></a>00104 <span class="preprocessor"></span> <span class="keywordtype">float</span> r=ip[0]; <span class="keywordtype">float</span> s=ip[1]; <span class="keywordtype">float</span> t=ip[2];
<a name="l00105"></a>00105 <span class="comment">//FIXME temporarily using hexahedron shape f for speed</span>
<a name="l00106"></a>00106 shapef[0] = (1 - r) * (1 - s) * (1 - t) * 0.125;
<a name="l00107"></a>00107 shapef[1] = (1 + r) * (1 - s) * (1 - t) * 0.125;
<a name="l00108"></a>00108 shapef[2] = (1 + r) * (1 + s) * (1 - t) * 0.125;
<a name="l00109"></a>00109 shapef[3] = (1 - r) * (1 + s) * (1 - t) * 0.125;
<a name="l00110"></a>00110 shapef[4] = (1 - r) * (1 - s) * (1 + t) * 0.125;
<a name="l00111"></a>00111 shapef[5] = (1 + r) * (1 - s) * (1 + t) * 0.125;
<a name="l00112"></a>00112 shapef[6] = (1 + r) * (1 + s) * (1 + t) * 0.125;
<a name="l00113"></a>00113 shapef[7] = (1 - r) * (1 + s) * (1 + t) * 0.125;
<a name="l00114"></a>00114 <span class="preprocessor">#else</span>
<a name="l00115"></a>00115 <span class="preprocessor"></span> shapef = element-><a class="code" href="classitk_1_1fem_1_1Element.html#a7b70476cd08c4fd56852599616b4675">ShapeFunctions</a>(ip);
<a name="l00116"></a>00116 <span class="preprocessor">#endif</span>
<a name="l00117"></a>00117 <span class="preprocessor"></span> }
<a name="l00118"></a>00118 <span class="keywordflow">else</span> <span class="keywordflow">if</span> (ImageDimension==2)
<a name="l00119"></a>00119 {
<a name="l00120"></a>00120 shapef = element-><a class="code" href="classitk_1_1fem_1_1Element.html#a7b70476cd08c4fd56852599616b4675">ShapeFunctions</a>(ip);
<a name="l00121"></a>00121 }
<a name="l00122"></a>00122 <span class="keywordtype">float</span> solval,posval;
<a name="l00123"></a>00123 detJ=element-><a class="code" href="classitk_1_1fem_1_1Element.html#6486f1c955e238628471d669d0c75357">JacobianDeterminant</a>(ip);
<a name="l00124"></a>00124
<a name="l00125"></a>00125 <span class="keywordflow">for</span>(<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> f=0; f<ImageDimension; f++)
<a name="l00126"></a>00126 {
<a name="l00127"></a>00127 solval=0.0;
<a name="l00128"></a>00128 posval=0.0;
<a name="l00129"></a>00129 <span class="keywordflow">for</span>(<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n=0; n<Nnodes; n++)
<a name="l00130"></a>00130 {
<a name="l00131"></a>00131 posval += shapef[n]*((element-><a class="code" href="classitk_1_1fem_1_1Element.html#cd805333750c8f7aaf73443b7831270f">GetNodeCoordinates</a>(n))[f]);
<a name="l00132"></a>00132 solval += shapef[n] * S-><a class="code" href="classitk_1_1fem_1_1Solution.html#4d4ea59c159d3d809bd5f7cb83f44818">GetSolutionValue</a>( element-><a class="code" href="classitk_1_1fem_1_1Element.html#3df768beea4e1d0c669fc7203350a8ba">GetNode</a>(n)->GetDegreeOfFreedom(f) , TotalSolutionIndex);
<a name="l00133"></a>00133 }
<a name="l00134"></a>00134 gsol[f]=solval;
<a name="l00135"></a>00135 gip[f]=posval;
<a name="l00136"></a>00136 }
<a name="l00137"></a>00137
<a name="l00138"></a>00138 <span class="comment">// Adjust the size of a force vector returned from the load object so</span>
<a name="l00139"></a>00139 <span class="comment">// that it is equal to the number of DOFs per node. If the Fg returned</span>
<a name="l00140"></a>00140 <span class="comment">// a vector with less dimensions, we add zero elements. If the Fg</span>
<a name="l00141"></a>00141 <span class="comment">// returned a vector with more dimensions, we remove the extra dimensions.</span>
<a name="l00142"></a>00142 force.fill(0.0);
<a name="l00143"></a>00143
<a name="l00144"></a>00144 force=l0->Fe(gip,gsol);
<a name="l00145"></a>00145 <span class="comment">// Calculate the equivalent nodal loads</span>
<a name="l00146"></a>00146 <span class="keywordflow">for</span>(<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> n=0; n<Nnodes; n++)
<a name="l00147"></a>00147 {
<a name="l00148"></a>00148 <span class="keywordflow">for</span>(<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> d=0; d<Ndofs; d++)
<a name="l00149"></a>00149 {
<a name="l00150"></a>00150 <a class="code" href="classitk_1_1fem_1_1Element.html#c55f879c15daa302723012b7f9d21b94">itk::fem::Element::Float</a> temp=shapef[n]*force[d]*w*detJ;
<a name="l00151"></a>00151 Fe[n*Ndofs+d] += temp;
<a name="l00152"></a>00152 }
<a name="l00153"></a>00153 }
<a name="l00154"></a>00154
<a name="l00155"></a>00155 }
<a name="l00156"></a>00156
<a name="l00157"></a>00157 }
<a name="l00158"></a>00158
<a name="l00159"></a>00159 };
<a name="l00160"></a>00160
<a name="l00161"></a>00161
<a name="l00162"></a>00162 <span class="keyword">template</span><<span class="keyword">class</span> TLoadClass>
<a name="l00163"></a>00163 <span class="keyword">const</span> <span class="keywordtype">bool</span> ImageMetricLoadImplementation<TLoadClass>::m_Registered = <span class="keyword">false</span>;
<a name="l00164"></a>00164
<a name="l00165"></a>00165
<a name="l00166"></a>00166 }} <span class="comment">// end namespace itk::fem</span>
<a name="l00167"></a>00167
<a name="l00168"></a>00168 <span class="preprocessor">#endif // #ifndef __itkFEMImageMetricLoadImplementation_h</span>
<a name="l00169"></a>00169 <span class="preprocessor"></span>
</pre></div></div>
<hr><address><small>
Generated at Thu May 7 22:41:57 2009 for ITK by <a href="http://www.stack.nl/~dimitri/doxygen/index.html"> <img
src="http://www.stack.nl/~dimitri/doxygen/doxygen.png" alt="doxygen"
align="middle" border=0 width=110 height=53>
</a> 1.5.9 written by <a href="mailto:dimitri@stack.nl">Dimitri van Heesch</a>,
© 1997-2000</small></address>
</body>
</html>
distrib
>
Mandriva
>
2010.0
>
i586
>
media
>
contrib-release
>
by-pkgid
>
d5e74628f0e673bb8680aebce32b2c04
>
files
>
15713
