itkfemelementbase.cxx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkFEMElementBase.cxx,v $
  Language:  C++
  Date:      $Date: 2007-12-31 18:35:16 $
  Version:   $Revision: 1.38 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/

// disable debug warnings in MS compiler
#ifdef _MSC_VER
#pragma warning(disable: 4786)
#endif 

#include "itkFEMElementBase.h"
#include "vnl/algo/vnl_svd.h"
#include "vnl/algo/vnl_qr.h"   



namespace itk {
namespace fem {




#ifdef FEM_BUILD_VISUALIZATION

/** Global scale factor for drawing on the DC */
double Element::DC_Scale=1000.0;

/** Global scale factor for drawing on the DC */
double &Element::Node::DC_Scale=Element::DC_Scale;


/**
 * draws the node on DC
 */
void Element::Node::Draw(CDC* pDC, Solution::ConstPointer sol) const 
{
  // We can only draw 2D nodes here
  if(m_coordinates.size()!=2) { return; }

  
  // Normally we draw a white circle.
  CPen pen(PS_SOLID, 0, (COLORREF) RGB(0,0,0) );
  CBrush brush( RGB(255,255,255) );

  CPen* pOldpen=pDC->SelectObject(&pen);
  CBrush* pOldbrush=pDC->SelectObject(&brush);

  int x1=m_coordinates[0]*DC_Scale;
  int y1=m_coordinates[1]*DC_Scale;
  x1+=sol->GetSolutionValue(this->GetDegreeOfFreedom(0))*DC_Scale;
  y1+=sol->GetSolutionValue(this->GetDegreeOfFreedom(1))*DC_Scale;

  CPoint r1=CPoint(0,0);
  CPoint r=CPoint(5,5);

  pDC->DPtoLP(&r1);
  pDC->DPtoLP(&r);
  r=r-r1;

  pDC->Ellipse(x1-r.x, y1-r.y, x1+r.x, y1+r.y);

  pDC->SelectObject(pOldbrush);
  pDC->SelectObject(pOldpen);

}

#endif




/*
 * Read the Node from the input stream
 */
void Element::Node::Read(  std::istream& f, void* info )
{
  unsigned int n;

  /*
   * First call the parent's read function
   */
  Superclass::Read(f,info);

  /*
   * Read and set node coordinates
   */
  this->SkipWhiteSpace(f); f>>n; if(!f) goto out;
  this->m_coordinates.set_size(n);
  this->SkipWhiteSpace(f); f>>this->m_coordinates; if(!f) goto out;

out:

  if( !f )
  {
    throw FEMExceptionIO(__FILE__,__LINE__,"Element::Node::Read()","Error reading FEM node!");
  }

}




/*
 * Write the Node to the output stream
 */
void Element::Node::Write( std::ostream& f ) const 
{
  /**
   * First call the parent's write function
   */
  Superclass::Write(f);

  /**
   * Write actual data (node, and properties numbers)
   */
  
  /* write the value of dof */
  f<<"\t"<<this->m_coordinates.size();
  f<<" "<<this->m_coordinates<<"\t% Node coordinates"<<"\n";

  /** check for errors */
  if (!f)
  {
    throw FEMExceptionIO(__FILE__,__LINE__,"Element::Node::Write()","Error writing FEM node!");
  }
}




//////////////////////////////////////////////////////////////////////////
/*
 * Physics of a problem.
 */

/*
 * Compute and return element stiffnes matrix (Ke) in global coordinate
 * system.
 * The base class provides a general implementation which only computes
 *
 *     b   T
 * int    B(x) D B(x) dx
 *     a
 *
 * using the Gaussian numeric integration method. The function calls
 * GetIntegrationPointAndWeight() / GetNumberOfIntegrationPoints() to obtain
 * the integration points. It also calls the GetStrainDisplacementMatrix()
 * and GetMaterialMatrix() member functions.
 *
 * \param Ke Reference to the resulting stiffnes matrix.
 *
 * \note This is a very generic implementation of the stiffness matrix
 *       that is suitable for any problem/element definition. A specifc
 *       element may override this implementation with its own simple one.
 */
void Element::GetStiffnessMatrix(MatrixType& Ke) const
{
  // B and D matrices
  MatrixType B,D;
  this->GetMaterialMatrix( D );

  unsigned int Nip=this->GetNumberOfIntegrationPoints();

  VectorType ip;
  Float w;
  MatrixType J;
  MatrixType shapeDgl;
  MatrixType shapeD;

  // Calculate the contribution of 1st int. point to initialize
  // the Ke matrix to proper number of elements.
  this->GetIntegrationPointAndWeight(0,ip,w);
  this->ShapeFunctionDerivatives(ip,shapeD);
  this->Jacobian(ip,J,&shapeD);
  this->ShapeFunctionGlobalDerivatives(ip,shapeDgl,&J,&shapeD);

  this->GetStrainDisplacementMatrix( B, shapeDgl );
  Float detJ=this->JacobianDeterminant( ip, &J );
  Ke=detJ*w*B.transpose()*D*B; // FIXME: write a more efficient way of computing this.

  // Add contributions of other int. points to the Ke
  for(unsigned int i=1; i<Nip; i++)
  {
    this->GetIntegrationPointAndWeight(i,ip,w);
    this->ShapeFunctionDerivatives(ip,shapeD);
    this->Jacobian(ip,J,&shapeD);
    this->ShapeFunctionGlobalDerivatives(ip,shapeDgl,&J,&shapeD);

    this->GetStrainDisplacementMatrix( B, shapeDgl );
    detJ=this->JacobianDeterminant( ip, &J );
    Ke+=detJ*w*B.transpose()*D*B; // FIXME: write a more efficient way of computing this.
  }
}

Element::VectorType Element::GetStrainsAtPoint(const VectorType& pt, const Solution& sol, unsigned int index) const
  // NOTE: pt should be in local coordinates already
{
  MatrixType B;
  VectorType e, u;
  MatrixType J, shapeD, shapeDgl;
  
  this->ShapeFunctionDerivatives(pt, shapeD);
  this->Jacobian(pt, J, &shapeD);
  this->ShapeFunctionGlobalDerivatives(pt, shapeDgl, &J, &shapeD);
  this->GetStrainDisplacementMatrix(B, shapeDgl);
  
  u = this->InterpolateSolution(pt, sol, index);
  
  e = B*u;

  return e;
}

Element::VectorType Element::GetStressesAtPoint(const VectorType& itkNotUsed(pt),
                                                const VectorType& e,
                                                const Solution& itkNotUsed(sol),
                                                unsigned int itkNotUsed(index)) const
  // NOTE: pt should be in local coordinates already
{
  MatrixType D;
  VectorType sigma;

  this->GetMaterialMatrix(D);
  
  sigma = D*e;

  return sigma;
}

void Element::GetLandmarkContributionMatrix(float eta, MatrixType& Le) const
{
  // Provides the contribution of a landmark to the element stiffness
  // matrix
  Le = MatrixType( this->GetNumberOfDegreesOfFreedom(), this->GetNumberOfDegreesOfFreedom(), 0.0 );

  const unsigned int NnDOF=this->GetNumberOfDegreesOfFreedomPerNode();
  const unsigned int Nnodes=this->GetNumberOfNodes();
  const unsigned int NDOF = GetNumberOfDegreesOfFreedom();
  const unsigned int Nip=this->GetNumberOfIntegrationPoints(0);

  Le.set_size(NDOF,NDOF); // resize the target matrix object
  Le.fill(0.0);

  Float w;
  VectorType ip, shape;

  for(unsigned int i=0; i<Nip; i++)
  {
    this->GetIntegrationPointAndWeight(i,ip,w,0);
    shape=this->ShapeFunctions(ip);
    
    for(unsigned int ni=0; ni<Nnodes; ni++)
    {
      for(unsigned int nj=0; nj<Nnodes; nj++)
      {
        Float m=w*shape[ni]*shape[nj];
        for(unsigned int d=0; d<NnDOF; d++)
        {
          Le[ni*NnDOF+d][nj*NnDOF+d]+=m;
        }
      }
    }
  }

  Le = Le / (eta );
}


Element::Float Element::GetElementDeformationEnergy( MatrixType& LocalSolution ) const
{

  MatrixType U;

  Element::MatrixType Ke;
  this->GetStiffnessMatrix(Ke);

  U=LocalSolution.transpose()*Ke*LocalSolution;

  return U[0][0];