itkfemlinearsystemwrapper.h

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

   * Virtual function to add a value of specific element of the solution
   * vector.
   * \param i element Index in solution vector
   * \param value new value of the element
   * \param solutionIndex index of solution vector to add value to
   */
  virtual void AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex = 0) = 0;

  /**
   * Solves the linear system and creates the solution vector, which can later
   * be accessed via GetSolutionValue(i,SolutionIndex) member function. Here all the major processing is
   * done with calls to external numeric library.
   * \note This function can only be called after the linear system was
   *       properly assembled.
   */
  virtual void Solve(void) = 0;

  /** 
   * Swaps access indices of any 2 matrices in the linear system
   * \param matrixIndex1 index of a matrix to swap
   * \param matrixIndex2 index of matrix to swap with
   */
  virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2) = 0;

  /** 
   * Copies the content of source matrix to destination matrix. Any existing
   * data in destination matrix is overwritten.
   *
   * \param matrixIndex1 index of a matrix that will be copied
   * \param matrixIndex2 index of matrix to copy to
   */
  virtual void CopyMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);

  /** 
   * Swaps access indices of any 2 vectors in the linear system
   * \param vectorIndex1 index of a vector to swap
   * \param vectorIndex2 index of vector to swap with
   */
  virtual void SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2) = 0;

  /** 
   * Swaps access indices of any 2 solution vectors in the linear system
   * \param solutionIndex1 index of a solution vector to swap
   * \param solutionIndex2 index of solution vector to swap with
   */
  virtual void SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2) = 0;


  /**
   * Multiplies all elements of a matrix by a scalar
   * \param scale scalar to multiply all matrix values by
   * \param matrixIndex index of matrix to modify
   */
  virtual void ScaleMatrix(Float scale, unsigned int matrixIndex = 0);


  /**
   * Multiplies all elements of a vector by a scalar
   * \param scale scalar to multiply all vector values by
   * \param vectorIndex index of vector to modify
   */
  void ScaleVector(Float scale, unsigned int vectorIndex = 0);


  /**
   * Multiplies all elements of a solution by a scalar
   * \param scale scalar to multiply all solution values by
   * \param solutionIndex index of solution to modify
   */
  void ScaleSolution(Float scale, unsigned int solutionIndex = 0);

  /**
   * Perform a matrix*matrix operation and store the result in the linear system
   * \param leftMatrixIndex index of left matrix
   * \param rightMatrixIndex index of right matrix
   * \param resultMatrixIndex index of matrix where solution is stored
   */
  virtual void MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex) = 0;

  /** 
   * Adds two matrices storing the result in the first matrix.
   *
   * \param matrixIndex1 index of a matrix to add the other matrix to
   * \param matrixIndex2 index of matrix to add
   */
  virtual void AddMatrixMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);

  /** 
   * Adds two vectors storing the result in the first vector.
   *
   * \param vectorIndex1 index of a vector to add the other vector to
   * \param vectorIndex2 index of vector to add
   */
  virtual void AddVectorVector(unsigned int vectorIndex1, unsigned int vectorIndex2);

  /**
   * Perform a matrix*vector operation and store the result in the linear system
   * \param matrixIndex index of matrix to multiply
   * \param vectorIndex index of vector to multiply
   * \param resultVectorIndex index of vector where result is store
   */
  virtual void MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);

  /**
   * Copy a solution vector to a vector
   * \param solutionIndex index of solution vector to copy
   * \param vectorIndex index of vector to copy solution to
   */
  virtual void CopySolution2Vector(unsigned int solutionIndex, unsigned int vectorIndex) = 0;

  /**
   * Copy a vector to a solution vector
   * \param vectorIndex index of a vector to copy
   * \param solutionIndex index of a solution to copy the solution to
   */
  virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex) = 0;
  /**
   * Copy a vector
   * \param vectorSource index of a vector to copy
   * \param vectorDestination index to copy the vector to
   */
  virtual void CopyVector(unsigned int vectorSource, unsigned int vectorDestination);

  /**
   * Remove all zeros from a matrix 
   * \param matrixIndex index of matrix to remove zeros from
   * \param tempMatrixIndex index of matrix to use for temp storage space
   * \note an extra matrix must be allocated by the solver in order to use this method
   */
  virtual void OptimizeMatrixStorage(unsigned int matrixIndex, unsigned int tempMatrixIndex);

  /**
   * Reorder the Degrees of Freedom in order to reduce bandwidth of matrix
   * \param matrixIndex index of matrix to examine
   * \param newNumbering vector of new degree of freedom ordering
   */
  virtual void ReverseCuthillMckeeOrdering(ColumnArray& newNumbering, unsigned int matrixIndex = 0);

  /*
   * Sets the function used to prepare the primary system matrix for numerical solving
   * \param SetupFunction pointer to function that stores the matrix to 
   * solve in the 0 matrix of the linear system
   */
  /*
  void SetPrimaryMatrixSetupFunction(void (*SetupFunction)(LinearSystemWrapper *lsw))
  { 
    m_PrimaryMatrixSetupFunction = SetupFunction; 
  }
  */


  /*
   * Sets the function used to prepare the primary system vector for numerical solving
   * \param SetupFunction pointer to function that stores the vector to
   * solve in the 0 vector of the linear system
   */
  /*
  void SetPrimaryVectorSetupFunction(void (*SetupFunction)(LinearSystemWrapper *lsw))
  { 
    m_PrimaryVectorSetupFunction = SetupFunction; 
  }
  */

  /*
   * Sets the function used to prepare the primary system solution for numerical solving
   * \param SetupFunction pointer to function that stores the solution 
   * in the 0 solution vector of the linear system
   */
  /*
  void SetPrimarySolutionSetupFunction(void (*SetupFunction)(LinearSystemWrapper *lsw))
  { 
    m_PrimarySolutionSetupFunction = SetupFunction; 
  }
  */

protected:

  /** Order of linear system */
  unsigned int m_Order;

  /**
   * Number of matrices used by system 
   */
  unsigned int m_NumberOfMatrices;

  /**
   * Number of vectors used by system 
   */
  unsigned int m_NumberOfVectors;

  /**
   * Number of solutions used by system 
   */
  unsigned int m_NumberOfSolutions;

  /*
   * Function used to prepare primary matrix for numerical solving 
   */
  //void (*m_PrimaryMatrixSetupFunction)(LinearSystemWrapper *lsw);

  /*
   * Function used to prepare primary vector for numerical solving 
   */
  /* void (*m_PrimaryVectorSetupFunction)(LinearSystemWrapper *lsw);*/

  /*
   * Function used to prepare primary matrix for numerical solving 
   */
  /* void (*m_PrimarySolutionSetupFunction)(LinearSystemWrapper *lsw); */

private:

  /**
   * matrix reordering utility
   */
  void CuthillMckeeOrdering(ColumnArray& newNumbering, int startingRow, unsigned int matrixIndex = 0);

  void FollowConnectionsCuthillMckeeOrdering(unsigned int rowNumber, ColumnArray& rowDegree, ColumnArray& newNumbering, unsigned int nextRowNumber, unsigned int matrixIndex = 0);

  /** Copy constructor is not allowed. */
  LinearSystemWrapper(const LinearSystemWrapper&);

  /** Asignment operator is not allowed. */
  const LinearSystemWrapper& operator= (const LinearSystemWrapper&);

};



class FEMExceptionLinearSystem : public FEMException
{
public:
  /**
   * Constructor. In order to construct this exception object, four parameters
   * must be provided: file, lineNumber, location and a detailed description
   * of the exception.
   */
  FEMExceptionLinearSystem(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription);
 
  /** Virtual destructor needed for subclasses. Has to have empty throw(). */
  virtual ~FEMExceptionLinearSystem() throw() {}
  
  /** Type related information. */
  itkTypeMacro(FEMExceptionLinearSystem,FEMException);
  
};

class FEMExceptionLinearSystemBounds : public FEMException
{
public:
  /**
   * Constructor. In order to construct this exception object, five parameters
   * must be provided: file, lineNumber, location and a detailed description
   * of the exception, and the invalid index
   */
  FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription, unsigned int index1);
 
  /**
   * Constructor. In order to construct this exception object, six parameters
   * must be provided: file, lineNumber, location and a detailed description
   * of the exception, the first index, and the second index   */
  FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription, unsigned int index1, unsigned int index2);

  /** Virtual destructor needed for subclasses. Has to have empty throw(). */
  virtual ~FEMExceptionLinearSystemBounds() throw() {}
  
  /** Type related information. */
  itkTypeMacro(FEMExceptionLinearSystem,FEMException);
  
};

}} // end namespace itk::fem

#endif // #ifndef __itkFEMLinearSystemWrapper_h