itkfemsolver.cxx

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

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

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

  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 "itkFEMSolver.h"

#include "itkFEMLoadNode.h"
#include "itkFEMLoadElementBase.h"
#include "itkFEMElementBase.h"
#include "itkFEMLoadBC.h"
#include "itkFEMLoadBCMFC.h"
#include "itkFEMLoadLandmark.h"

#include "itkImageRegionIterator.h"

#include <algorithm>

namespace itk {
namespace fem {

/*
 * Default constructor for Solver class
 */
Solver::Solver() : NGFN(0), NMFC(0)
{
  this->SetLinearSystemWrapper(&m_lsVNL);
}

void Solver::Clear( void )
{
  this->el.clear();
  this->node.clear();
  this->mat.clear();
  this->load.clear();

  this->NGFN=0;
  this->NMFC=0;
  this->SetLinearSystemWrapper(&m_lsVNL);
}


/*
 * Change the LinearSystemWrapper object used to solve
 * system of equations.
 */
void Solver::SetLinearSystemWrapper(LinearSystemWrapper::Pointer ls)
{
  m_ls=ls; // update the pointer to LinearSystemWrapper object

  this->InitializeLinearSystemWrapper();
}


void Solver::InitializeLinearSystemWrapper(void)
{
  // set the maximum number of matrices and vectors that
  // we will need to store inside.
  m_ls->SetNumberOfMatrices(1);
  m_ls->SetNumberOfVectors(2);
  m_ls->SetNumberOfSolutions(1);
}




/*
 * Reads the whole system (nodes, materials and elements) from input stream
 */
void Solver::Read(std::istream& f) {

  // clear all arrays - uncomment these 4 lines if you want to use the
  // Windows-based visualization
  //   el.clear();
  //   node.clear();
  //   mat.clear();
  //   load.clear();

  // Initialize the pointers to arrays in ReadInfoType object to the
  // arrays in solver object.
  ReadInfoType info(&this->node,&this->el,&this->mat);

  FEMLightObject::Pointer o;
  /* then we start reading objects from stream */
  do
    {
    o=FEMLightObject::CreateFromStream(f,&info);
    /*
     * If CreateFromStream returned 0, we're ok. That was the signal
     * for the end of stream. Just continue reading... and consequently
     * exit the do loop.
     */
    if (!o) { continue; }

    /*
     * Find out what kind of object did we read from stream
     * and store it in the appropriate array
     */
    if ( Node::Pointer o1=dynamic_cast<Node*>(&*o) )
      {
      node.push_back(FEMP<Node>(o1));
      continue;
      }
    if ( Material::Pointer o1=dynamic_cast<Material*>(&*o) )
      {
      mat.push_back(FEMP<Material>(o1));
      continue;
      }
    if ( Element::Pointer o1=dynamic_cast<Element*>(&*o) )
      {
      el.push_back(FEMP<Element>(o1));
      continue;
      }
    if ( Load::Pointer o1=dynamic_cast<Load*>(&*o) )
      {
      load.push_back(FEMP<Load>(o1));
      continue;
      }

    /*
     * If we got here, something strange was in the file...
     */

    // first we delete the allocated object
#ifndef FEM_USE_SMART_POINTERS
    delete o;
#endif
    o=0;

    // then we throw an exception
    throw FEMExceptionIO(__FILE__,__LINE__,"Solver::Read()","Error reading FEM problem stream!");

    } while ( o );

}


/*
 * Writes everything (nodes, materials and elements) to output stream
 */
void Solver::Write( std::ostream& f ) {

  for(NodeArray::iterator i=node.begin(); i!=node.end(); i++)
    {
    (*i)->Write(f);
    }
  f<<"\n<END>  % End of nodes\n\n";

  for(MaterialArray::iterator i=mat.begin(); i!=mat.end(); i++)
    {
    (*i)->Write(f);
    }
  f<<"\n<END>  % End of materials\n\n";

  for(ElementArray::iterator i=el.begin(); i!=el.end(); i++)
    {
    (*i)->Write(f);
    }
  f<<"\n<END>  % End of elements\n\n";

  for(LoadArray::iterator i=load.begin(); i!=load.end(); i++)
    {
    (*i)->Write(f);
    }
  f<<"\n<END>  % End of loads\n\n";

}




/*
 * Assign a global freedom number to each DOF in a system.
 */
void Solver::GenerateGFN() {

  // Clear the list of elements and global freedom numbers in nodes
  // FIXME: should be removed once Mesh is there
  for(NodeArray::iterator n=node.begin(); n!=node.end(); n++)
    {
    (*n)->m_elements.clear();
    (*n)->ClearDegreesOfFreedom();
    }

  for(ElementArray::iterator e=el.begin(); e!=el.end(); e++) // step over all elements
    {

    // Add the elemens in the nodes list of elements
    // FIXME: should be removed once Mesh is there
    unsigned int Npts=(*e)->GetNumberOfNodes();
    for(unsigned int pt=0; pt<Npts; pt++)
      {
      (*e)->GetNode(pt)->m_elements.insert(*e);
      }
    }


  /*
   * Assign new ID to every DOF in a system
   */

  // Start numbering DOFs from 0
  NGFN=0;

  // Step over all elements
  for(ElementArray::iterator e=el.begin(); e!=el.end(); e++)
    {
    // FIXME: Write a code that checks if two elements are compatible, when they share a node
    for(unsigned int n=0; n<(*e)->GetNumberOfNodes(); n++)
      {
      for(unsigned int dof=0; dof<(*e)->GetNumberOfDegreesOfFreedomPerNode(); dof++)
        {
        if( (*e)->GetNode(n)->GetDegreeOfFreedom(dof)==Element::InvalidDegreeOfFreedomID )
          {
          (*e)->GetNode(n)->SetDegreeOfFreedom(dof,NGFN);
          NGFN++;
          }
        }
      }
    } // end for e

//  NGFN=Element::GetGlobalDOFCounter()+1;
  if (NGFN>0) return;  // if we got 0 DOF, somebody forgot to define the system...

}




/*
 * Assemble the master stiffness matrix (also apply the MFCs to K)
 */
void Solver::AssembleK()
{

  // if no DOFs exist in a system, we have nothing to do
  if (NGFN<=0) return;

  NMFC=0;  // reset number of MFC in a system

  /*
   * Before we can start the assembly procedure, we need to know,
   * how many boundary conditions if form of MFCs are there in a system.
   */

  // search for MFC's in Loads array, because they affect the master stiffness matrix
  for(LoadArray::iterator l=load.begin(); l!=load.end(); l++)
    {
    if ( LoadBCMFC::Pointer l1=dynamic_cast<LoadBCMFC*>( &(*(*l))) ) {
    // store the index of an LoadBCMFC object for later
    l1->Index=NMFC;
    // increase the number of MFC
    NMFC++;
    }
    }

  /*
   * Now we can assemble the master stiffness matrix from
   * element stiffness matrices.
   *
   * Since we're using the Lagrange multiplier method to apply the MFC,
   * each constraint adds a new global DOF.
   */
  this->InitializeMatrixForAssembly(NGFN+NMFC);

  /*
   * Step over all elements
   */
  for(ElementArray::iterator e=el.begin(); e!=el.end(); e++)
    {
    // Call the function that actually moves the element matrix
    // to the master matrix.
    this->AssembleElementMatrix(&**e);
    }

  /*
   * Step over all the loads again to add the landmark contributions
   * to the appropriate place in the stiffness matrix
   */
  for(LoadArray::iterator l2=load.begin(); l2!=load.end(); l2++)
    {
    if ( LoadLandmark::Pointer l3=dynamic_cast<LoadLandmark*>( &(*(*l2))) )
      {
      l3->AssignToElement(&el);
      Element::Pointer ep = const_cast<Element*>( l3->el[0] );
      this->AssembleLandmarkContribution( ep , l3->eta );
      }
    }

  this->FinalizeMatrixAfterAssembly();

}




void Solver::InitializeMatrixForAssembly(unsigned int N)
{
  // We use LinearSystemWrapper object, to store the K matrix.
  this->m_ls->SetSystemOrder(N);
  this->m_ls->InitializeMatrix();
}


void Solver::AssembleLandmarkContribution(Element::Pointer e, float eta)
{
  // Copy the element "landmark" matrix for faster access.
  Element::MatrixType Le;
  e->GetLandmarkContributionMatrix(eta, Le);

  // ... same for number of DOF
  int Ne=e->GetNumberOfDegreesOfFreedom();

  // step over all rows in element matrix
  for(int j=0; j<Ne; j++)
    {
    // step over all columns in element matrix
    for(int k=0; k<Ne; k++)
      {
      // error checking. all GFN should be =>0 and <NGFN
      if ( e->GetDegreeOfFreedom(j) >= NGFN ||
           e->GetDegreeOfFreedom(k) >= NGFN  )
        {
        throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::AssembleLandmarkContribution()","Illegal GFN!");
        }

      /*
       * Here we finaly update the corresponding element
       * in the master stiffness matrix. We first check if
       * element in Le is zero, to prevent zeros from being
       * allocated in sparse matrix.
       */
      if ( Le[j][k]!=Float(0.0) )
        {
        this->m_ls->AddMatrixValue( e->GetDegreeOfFreedom(j), e->GetDegreeOfFreedom(k), Le[j][k] );
        }
      }
    }
}



void Solver::AssembleElementMatrix(Element::Pointer e)
{
  // Copy the element stiffness matrix for faster access.
  Element::MatrixType Ke;
  e->GetStiffnessMatrix(Ke);

  // ... same for number of DOF
  int Ne=e->GetNumberOfDegreesOfFreedom();

  // step over all rows in element matrix
  for(int j=0; j<Ne; j++)
    {
    // step over all columns in element matrix
    for(int k=0; k<Ne; k++)
      {
      // error checking. all GFN should be =>0 and <NGFN
      if ( e->GetDegreeOfFreedom(j) >= NGFN ||
           e->GetDegreeOfFreedom(k) >= NGFN  )
        {
        throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::AssembleElementMatrix()","Illegal GFN!");
        }

      /*
       * Here we finaly update the corresponding element
       * in the master stiffness matrix. We first check if
       * element in Ke is zero, to prevent zeros from being
       * allocated in sparse matrix.
       */
      if ( Ke[j][k]!=Float(0.0) )
        {
        this->m_ls->AddMatrixValue( e->GetDegreeOfFreedom(j), e->GetDegreeOfFreedom(k), Ke[j][k] );
        }

      }

    }

}




/*
 * Assemble the master force vector
 */
void Solver::AssembleF(int dim) {

  // Vector that stores element nodal loads
  Element::VectorType Fe;

  // Type that stores IDs of fixed DOF together with the values to
  // which they were fixed.
  typedef std::map<Element::DegreeOfFreedomIDType,Float> BCTermType;
  BCTermType bcterm;

  /* if no DOFs exist in a system, we have nothing to do */
  if (NGFN<=0) return;

  /* Initialize the master force vector */
  m_ls->InitializeVector();

  /*
   * Convert the external loads to the nodal loads and
   * add them to the master force vector F.
   */
  for(LoadArray::iterator l=load.begin(); l!=load.end(); l++)
    {

    /*
     * Store a temporary pointer to load object for later,
     * so that we don't have to access it via the iterator
     */
    Load::Pointer l0=*l;

    /*
     * Pass the vector to the solution to the Load object.
     */
    l0->SetSolution(m_ls);

    /*
     * Here we only handle Nodal loads
     */
    if ( LoadNode::Pointer l1=dynamic_cast<LoadNode*>(&*l0) )
      {
      // yep, we have a nodal load
      
      // size of a force vector in load must match number of DOFs in node
      if ( (l1->F.size() % l1->m_element->GetNumberOfDegreesOfFreedomPerNode())!=0 )
        {
        throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::AssembleF()","Illegal size of a force vector in LoadNode object!");
        }

      // we simply copy the load to the force vector
      for(unsigned int d=0; d < (l1->m_element->GetNumberOfDegreesOfFreedomPerNode()); d++)
        {
        Element::DegreeOfFreedomIDType dof=l1->m_element->GetNode(l1->m_pt)->GetDegreeOfFreedom(d);
        // error checking
        if ( dof >= NGFN )
          {
          throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::AssembleF()","Illegal GFN!");
          }

        /*
         * If using the extra dim parameter, we can apply the force to
         * different isotropic dimension.
         *
         * FIXME: We assume that the implementation of force vector
         * inside the LoadNode class is correct for given number of