itkFEMElementBase.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

#ifndef __itkFEMElementBase_h
#define __itkFEMElementBase_h

#include "itkFEMLightObject.h"
#include "itkFEMPArray.h"
#include "itkFEMMaterialBase.h"
#include "itkFEMSolution.h"
#include "itkVectorContainer.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_vector.h"

#include <set>
#include <vector>

namespace itk
{
namespace fem
{
class Element : public FEMLightObject
{
public:
  typedef Element                  Self;
  typedef FEMLightObject           Superclass;
  typedef SmartPointer<Self>       Pointer;
  typedef SmartPointer<const Self> ConstPointer;

  itkTypeMacro(Element, FEMLightObject);

  typedef double        Float;
  typedef unsigned long ElementIdentifier;

  // FIXME - Remove FEMPArray Type and replace with VectorContainer version
  typedef FEMPArray<Element>                                   ArrayType;
  typedef VectorContainer<ElementIdentifier, Element::Pointer> ArrayType1;

  typedef vnl_matrix<Float> MatrixType;

  typedef vnl_vector<Float> VectorType;

  typedef FEMLightObject    LoadType;
  typedef LoadType::Pointer LoadPointer;

  typedef unsigned int DegreeOfFreedomIDType;

  enum { InvalidDegreeOfFreedomID = 0xffffffff };


  class Node : public FEMLightObject
  {
  public:

    typedef Node                     Self;
    typedef FEMLightObject           Superclass;
    typedef SmartPointer<Self>       Pointer;
    typedef SmartPointer<const Self> ConstPointer;

    // itkNewMacro(Self);
    static Pointer New(void)
      {
        Pointer smartPtr = ::itk::ObjectFactory<Self>::Create();

        if( smartPtr.IsNull() )
          {
          smartPtr = static_cast<Pointer>(new Self);
          }
        smartPtr->UnRegister();
        return smartPtr;
      }

    itkTypeMacro(Node, FEMLightObject);

    virtual::itk::LightObject::Pointer CreateAnother(void) const
      {
        ::itk::LightObject::Pointer smartPtr;
        Pointer copyPtr = Self::New();

        copyPtr->m_coordinates = this->m_coordinates;
        copyPtr->m_dof = this->m_dof;
        copyPtr->m_elements = this->m_elements;
        copyPtr->SetGlobalNumber( this->GetGlobalNumber() );

        smartPtr = static_cast<Pointer>(copyPtr);

        return smartPtr;
      }

    typedef double Float;

    typedef FEMPArray<Self> ArrayType;

    Node()
      {
      }

    ~Node()
      {
        this->ClearDegreesOfFreedom();
        this->m_elements.clear();
      }

    const VectorType & GetCoordinates(void) const
      {
        return m_coordinates;
      }

    void SetCoordinates(const VectorType & coords)
      {
        m_coordinates = coords;
      }

    DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int i) const
      {
        if( i >= m_dof.size() )
          {
          return InvalidDegreeOfFreedomID;
          }
        return m_dof[i];
      }

    void SetDegreeOfFreedom(unsigned int i, DegreeOfFreedomIDType dof) const
      {
        if( i >= m_dof.size() )
          {
          m_dof.resize(i + 1, InvalidDegreeOfFreedomID);
          }
        m_dof[i] = dof;
      }

    virtual void ClearDegreesOfFreedom(void) const
      {
        m_dof.clear();
      }

  public:
    typedef std::set<Element *> SetOfElements;
    mutable SetOfElements m_elements;
  protected:
    virtual void PrintSelf(std::ostream& os, Indent indent) const
      {
        Superclass::PrintSelf(os, indent);
        // os << indent << "DOF: " << this->m_dof << std::endl;
        // os << indent << "Coordinates: " << this->m_coordinates << std::endl;
        // os << indent << "Elements: " << this->m_elements << std::endl;
      }

  private:
    VectorType m_coordinates;

    mutable std::vector<DegreeOfFreedomIDType> m_dof;
  };  // end class Node

// ////////////////////////////////////////////////////////////////////////
/*
 * Methods related to the physics of the problem.
 */

  virtual VectorType GetStrainsAtPoint(const VectorType & pt, const Solution & sol, unsigned int index) const;

  virtual VectorType GetStressesAtPoint(const VectorType & pt, const VectorType & e, const Solution & sol,
                                        unsigned int index) const;

  virtual void GetStiffnessMatrix(MatrixType & Ke) const;

  virtual Float GetElementDeformationEnergy(MatrixType & LocalSolution) const;

  virtual void GetMassMatrix(MatrixType & Me) const;

  virtual void GetLandmarkContributionMatrix(float eta, MatrixType & Le) const;

  virtual void GetStrainDisplacementMatrix(MatrixType & B, const MatrixType & shapeDgl) const = 0;

  virtual void GetMaterialMatrix(MatrixType & D) const = 0;

  virtual VectorType InterpolateSolution(const VectorType & pt,
                                         const Solution & sol,
                                         unsigned int solutionIndex = 0) const;

  virtual Float InterpolateSolutionN(const VectorType & pt, const Solution & sol, unsigned int f,
                                     unsigned int solutionIndex = 0) const;

  DegreeOfFreedomIDType GetDegreeOfFreedom(unsigned int local_dof) const
    {
      if( local_dof > this->GetNumberOfDegreesOfFreedom() )
        {
        return InvalidDegreeOfFreedomID;
        }
      return this->GetNode(local_dof /
                           this->GetNumberOfDegreesOfFreedomPerNode() )
        ->GetDegreeOfFreedom(local_dof % this->GetNumberOfDegreesOfFreedomPerNode() );
    }

  virtual Material::ConstPointer GetMaterial(void) const
    {
      return 0;
    }

  virtual void SetMaterial(Material::ConstPointer)
    {
    }

  // ////////////////////////////////////////////////////////////////////////
  virtual void GetIntegrationPointAndWeight(unsigned int i,
                                            VectorType & pt,
                                            Float & w,
                                            unsigned int order = 0) const = 0;

  virtual unsigned int GetNumberOfIntegrationPoints(unsigned int order = 0) const = 0;

  enum { gaussMaxOrder = 10 };

  static const Float gaussPoint[gaussMaxOrder + 1][gaussMaxOrder];

  static const Float gaussWeight[gaussMaxOrder + 1][gaussMaxOrder];

// ////////////////////////////////////////////////////////////////////////
/*
 * Methods related to the geometry of an element
 */

  typedef Node::ConstPointer NodeIDType;

  virtual unsigned int GetNumberOfNodes(void) const = 0;

  virtual NodeIDType GetNode(unsigned int n) const = 0;

  virtual void SetNode(unsigned int n, NodeIDType node) = 0;
  virtual void SetNode(unsigned int n, Node::Pointer node)
    {
      this->SetNode(n,NodeIDType(node.GetPointer()));
    }

  virtual const VectorType & GetNodeCoordinates(unsigned int n) const = 0;

  virtual VectorType GetGlobalFromLocalCoordinates(const VectorType & pt) const;

  virtual bool GetLocalFromGlobalCoordinates(const VectorType & globalPt, VectorType & localPt) const = 0;

  virtual unsigned int GetNumberOfSpatialDimensions() const = 0;

  virtual VectorType ShapeFunctions(const VectorType & pt) const = 0;

  virtual void ShapeFunctionDerivatives(const VectorType & pt, MatrixType & shapeD) const = 0;

  virtual void ShapeFunctionGlobalDerivatives(const VectorType & pt, MatrixType & shapeDgl, const MatrixType *pJ = 0,
                                              const MatrixType *pshapeD = 0) const;

  virtual void Jacobian(const VectorType & pt, MatrixType & J, const MatrixType *pshapeD = 0) const;

  virtual Float JacobianDeterminant(const VectorType & pt, const MatrixType *pJ = 0) const;

  virtual void JacobianInverse(const VectorType & pt, MatrixType & invJ, const MatrixType *pJ = 0) const;

  virtual unsigned int GetNumberOfDegreesOfFreedom(void) const
    {
      return this->GetNumberOfNodes() * this->GetNumberOfDegreesOfFreedomPerNode();
    }

  virtual std::vector<std::vector<int> > GetEdgeIds(void) const
    {
      return this->m_EdgeIds;
    }

  virtual unsigned int GetNumberOfDegreesOfFreedomPerNode(void) const = 0;

  virtual void PopulateEdgeIds(void) = 0;

protected:

  // to store edge connectivity data
  std::vector<std::vector<int> > m_EdgeIds;

  virtual void PrintSelf(std::ostream& os, Indent indent) const;

};

}
}  // end namespace itk::fem

#endif // #ifndef __itkFEMElementBase_h