itkfemsolver.cxx

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

         * dimensions
         */
        m_ls->AddVectorValue(dof , l1->F[d+l1->m_element->GetNumberOfDegreesOfFreedomPerNode()*dim]);
        }

      // that's all there is to DOF loads, go to next load in an array
      continue;
      }


    /*
     * Element loads...
     */
    if ( LoadElement::Pointer l1=dynamic_cast<LoadElement*>(&*l0) )
      {

      if ( !(l1->el.empty()) )
        {
        /*
         * If array of element pointers is not empty,
         * we apply the load to all elements in that array
         */
        for(LoadElement::ElementPointersVectorType::const_iterator i=l1->el.begin(); i!=l1->el.end(); i++)
          {

          const Element* el0=(*i);
          // Call the Fe() function of the element that we are applying the load to.
          // We pass a pointer to the load object as a paramater and a reference to the nodal loads vector.
          el0->GetLoadVector(Element::LoadPointer(l1),Fe);
          unsigned int Ne=el0->GetNumberOfDegreesOfFreedom();    // ... element's number of DOF
          for(unsigned int j=0; j<Ne; j++)    // step over all DOF
            {
            // error checking
            if ( el0->GetDegreeOfFreedom(j) >= NGFN )
              {
              throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::AssembleF()","Illegal GFN!");
              }

            // update the master force vector (take care of the correct isotropic dimensions)
            m_ls->AddVectorValue(el0->GetDegreeOfFreedom(j) , Fe(j+dim*Ne));
            }
          }

        } 
      else
        {

        /*
         * If the list of element pointers in load object is empty,
         * we apply the load to all elements in a system.
         */
        for(ElementArray::iterator e=el.begin(); e!=el.end(); e++) // step over all elements in a system
          {
          (*e)->GetLoadVector(Element::LoadPointer(l1),Fe);  // ... element's force vector
          unsigned int Ne=(*e)->GetNumberOfDegreesOfFreedom();    // ... element's number of DOF

          for(unsigned int j=0; j<Ne; j++)  // step over all DOF
            {
            if ( (*e)->GetDegreeOfFreedom(j) >= NGFN )
              {
              throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::AssembleF()","Illegal GFN!");
              }

            // update the master force vector (take care of the correct isotropic dimensions)
            m_ls->AddVectorValue((*e)->GetDegreeOfFreedom(j) , Fe(j+dim*Ne));
            
            }
          }
        }

      // skip to next load in an array
      continue;
      }

    /*
     * Handle boundary conditions in form of MFC loads are handled next.
     */
    if ( LoadBCMFC::Pointer l1=dynamic_cast<LoadBCMFC*>(&*l0) )
      {
      m_ls->SetVectorValue(NGFN+l1->Index , l1->rhs[dim]);

      // skip to next load in an array
      continue;
      }

    /*
     * Handle essential boundary conditions.
     */
    if ( LoadBC::Pointer l1=dynamic_cast<LoadBC*>(&*l0) )
      {
      
      // Here we just store the values of fixed DOFs. We can't set it here, because
      // it may be changed by other loads that are applied later.
      bcterm[ l1->m_element->GetDegreeOfFreedom(l1->m_dof) ]=l1->m_value[dim];
      
      // skip to the next load in an array
      continue;
      }
    
    /*
     * If we got here, we were unable to handle that class of Load object.
     * We do nothing...
     */
    

    }  // for(LoadArray::iterator l ... )

  /*
   * Adjust the master force vector for essential boundary
   * conditions as required.
   */
  if ( m_ls->IsVectorInitialized(1) )
    {
    // Add the vector generated by ApplyBC to the solution vector
    const unsigned int totGFN=NGFN+NMFC;
    for( unsigned int i=0; i<totGFN; i++ )
      {
      m_ls->AddVectorValue(i,m_ls->GetVectorValue(i,1));
      }

    }

  // Set the fixed DOFs to proper values
  for( BCTermType::iterator q=bcterm.begin(); q!=bcterm.end(); q++)
    {
    m_ls->SetVectorValue(q->first,q->second);
    }

}

/*
 * Decompose matrix using svd, qr, whatever ... if needed
 */
void Solver::DecomposeK()
{
}


/*
 * Solve for the displacement vector u
 */
void Solver::Solve()
{
  // Check if master stiffness matrix and master force vector were
  // properly initialized.
  if(!m_ls->IsMatrixInitialized())
    {
    throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::Solve()","Master stiffness matrix was not initialized!");
    }
  if(!m_ls->IsVectorInitialized())
    {
    throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::Solve()","Master force vector was not initialized!");
    }

  // Solve the system of linear equations
  m_ls->InitializeSolution();
  m_ls->Solve();
}

/*
 * Copy solution vector u to the corresponding nodal values, which are
 * stored in node objects). This is standard post processing of the solution.
 */
void Solver::UpdateDisplacements()
{
}

Solver::Float Solver::GetDeformationEnergy(unsigned int SolutionIndex)
{
  float U=0.0;
  Element::MatrixType LocalSolution;

  for(ElementArray::iterator e=el.begin(); e!=el.end(); e++)
    {
    unsigned int Ne=(*e)->GetNumberOfDegreesOfFreedom();
    LocalSolution.set_size(Ne,1);
    // step over all DOFs of element
    for(unsigned int j=0; j<Ne; j++)
      {
      LocalSolution[j][0]=m_ls->GetSolutionValue((*e)->GetDegreeOfFreedom(j),SolutionIndex);
      }

    U+=(*e)->GetElementDeformationEnergy(LocalSolution);
    }
  return U;
}

/*
 * Apply the boundary conditions to the system.
 */
void Solver::ApplyBC(int dim, unsigned int matrix)
{

  // Vector with index 1 is used to store force correctios for BCs
  m_ls->DestroyVector(1);

  /* Step over all Loads */
  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;


    /*
     * Apply boundary conditions in form of MFC loads.
     *
     * We add the multi freedom constraints contribution to the master
     * stiffness matrix using the lagrange multipliers. Basically we only
     * change the last couple of rows and columns in K.
     */
    if ( LoadBCMFC::Pointer c=dynamic_cast<LoadBCMFC*>(&*l0) )
      {
      /* step over all DOFs in MFC */
      for(LoadBCMFC::LhsType::iterator q=c->lhs.begin();
          q!=c->lhs.end();
          q++) {

      /* obtain the GFN of DOF that is in the MFC */
      Element::DegreeOfFreedomIDType gfn=
        q->m_element->GetDegreeOfFreedom(q->dof);

      /* error checking. all GFN should be =>0 and <NGFN */
      if ( gfn>=NGFN )
        {
        throw FEMExceptionSolution(__FILE__,__LINE__,"Solver::ApplyBC()","Illegal GFN!");
        }

      /* set the proper values in matster stiffnes matrix */
      this->m_ls->SetMatrixValue(gfn, NGFN+c->Index, q->value, matrix);
      this->m_ls->SetMatrixValue(NGFN+c->Index, gfn, q->value, matrix);  // this is a symetric matrix...

      }

      // skip to next load in an array
      continue;
      }

    /*
     * Apply essential boundary conditions
     */
    if ( LoadBC::Pointer c=dynamic_cast<LoadBC*>(&*l0) )
      {

      Element::DegreeOfFreedomIDType fdof = c->m_element->GetDegreeOfFreedom(c->m_dof);
      Float fixedvalue=c->m_value[dim];


      // Copy the corresponding row of the matrix to the vector that will
      // be later added to the master force vector.
      // NOTE: We need to copy the whole row first, and then clear it. This
      //       is much more efficient when using sparse matrix storage, than
      //       copying and clearing in one loop.

      // Get the column indices of the nonzero elements in an array.
      LinearSystemWrapper::ColumnArray cols;
      m_ls->GetColumnsOfNonZeroMatrixElementsInRow(fdof, cols, matrix);

      // Force vector needs updating only if DOF was not fixed to 0.0.
      if( fixedvalue!=0.0 )
        {
        // Initialize the master force correction vector as required
        if ( !this->m_ls->IsVectorInitialized(1) )
          {
          this->m_ls->InitializeVector(1);
          }

        // Step over each nonzero matrix element in a row
        for(LinearSystemWrapper::ColumnArray::iterator cc=cols.begin(); cc!=cols.end(); cc++)
          {
          // Get value from the stiffness matrix
          Float d=this->m_ls->GetMatrixValue(fdof, *cc, matrix);

          // Store the appropriate value in bc correction vector (-K12*u2)
          //
          // See http://titan.colorado.edu/courses.d/IFEM.d/IFEM.Ch04.d/IFEM.Ch04.pdf
          // chapter 4.1.3 (Matrix Forms of DBC Application Methods) for more info.
          this->m_ls->AddVectorValue(*cc,-d*fixedvalue,1);
          }
        }


      // Clear that row and column in master matrix
      for(LinearSystemWrapper::ColumnArray::iterator cc=cols.begin(); cc!=cols.end(); cc++)
        {
        this->m_ls->SetMatrixValue(fdof,*cc, 0.0, matrix);
        this->m_ls->SetMatrixValue(*cc,fdof, 0.0, matrix); // this is a symetric matrix
        }
      this->m_ls->SetMatrixValue(fdof,fdof, 1.0, matrix); // Set the diagonal element to one


      // skip to next load in an array
      continue;

      }

    } // end for LoadArray::iterator l
}

/*
 * Initialize the interpolation grid
 */
void Solver::InitializeInterpolationGrid(const VectorType& size, const VectorType& bb1, const VectorType& bb2)
{
  // Discard any old image object an create a new one
  m_InterpolationGrid=InterpolationGridType::New();

  // Set the interpolation grid (image) size, origin and spacing
  // from the given vectors, so that physical point of v1 is (0,0,0) and
  // phisical point v2 is (size[0],size[1],size[2]).
  InterpolationGridType::SizeType image_size={{1,1,1}};
  for(unsigned int i=0;i<size.size();i++)
    {
    image_size[i] = static_cast<InterpolationGridType::SizeType::SizeValueType>( size[i] );
    }
  Float image_origin[MaxGridDimensions]={0.0,0.0,0.0};
  for(unsigned int i=0;i<size.size();i++)
    {
    image_origin[i]=bb1[i];
    }
  Float image_spacing[MaxGridDimensions]={1.0,1.0,1.0};
  for(unsigned int i=0;i<size.size();i++)
    {
    image_spacing[i]=(bb2[i]-bb1[i])/(image_size[i]-1);
    }

  // All regions are the same
  m_InterpolationGrid->SetRegions(image_size);
  m_InterpolationGrid->Allocate();

  // Set origin and spacing
  m_InterpolationGrid->SetOrigin(image_origin);
  m_InterpolationGrid->SetSpacing(image_spacing);

  // Initialize all pointers in interpolation grid image to 0
  m_InterpolationGrid->FillBuffer(0);

  VectorType v1,v2;

  // Fill the interpolation grid with proper pointers to elements
  for(ElementArray::iterator e=el.begin(); e!=el.end(); e++)
    {
    // Get square boundary box of an element
    v1=(*e)->GetNodeCoordinates(0); // lower left corner
    v2=v1; // upper right corner

    const unsigned int NumberOfDimensions=(*e)->GetNumberOfSpatialDimensions();

    for(unsigned int n=1; n < (*e)->GetNumberOfNodes(); n++ )
      {
      const VectorType& v=(*e)->GetNodeCoordinates(n);
      for(unsigned int d=0; d < NumberOfDimensions; d++ )
        {
        if( v[d] < v1[d] )
          {
          v1[d]=v[d];
          }
        if( v[d] > v2[d] )
          {
          v2[d]=v[d];
          }
        }
      }

    // Convert boundary box corner points into discrete image indexes.
    InterpolationGridType::IndexType vi1,vi2;

    Point<Float,MaxGridDimensions> vp1,vp2,pt;
    for(unsigned int i=0;i<MaxGridDimensions;i++)
      {
      if ( i < NumberOfDimensions )
        {
        vp1[i]=v1[i];
        vp2[i]=v2[i];
        }
      else
        {
        vp1[i]=0.0;
        vp2[i]=0.0;
        }
      }

    // Obtain the Index of BB corner and check whether it is within image.
    // If it is not, we ignore the entire element.
    if(!m_InterpolationGrid->TransformPhysicalPointToIndex(vp1,vi1)) continue;
    if(!m_InterpolationGrid->TransformPhysicalPointToIndex(vp2,vi2)) continue;

    InterpolationGridType::SizeType region_size;
    for( unsigned int i=0; i<MaxGridDimensions; i++ )
      {
      region_size[i] = vi2[i]-vi1[i]+1;
      }
    InterpolationGridType::RegionType region(vi1,region_size);

    // Initialize the iterator that will step over all grid points within
    // element boundary box.
    ImageRegionIterator<InterpolationGridType> iter(m_InterpolationGrid,region);

    //
    // Update the element pointers in the points defined within the region.
    //
    VectorType global_point(NumberOfDimensions); // Point in the image as a vector.
    VectorType local_point; // Same point in local element coordinate system

    // Step over all points within the region
    for(iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
      {
      // Note: Iteratior is guarantied to be within image, since the
      //       elements with BB outside are skipped before.
      m_InterpolationGrid->TransformIndexToPhysicalPoint(iter.GetIndex(),pt);
      for(unsigned int d=0; d<NumberOfDimensions; d++)
        {
        global_point[d]=pt[d];
        }

      // If the point is within the element, we update the pointer at
      // this point in the interpolation grid image.
      if( (*e)->GetLocalFromGlobalCoordinates(global_point,local_point) )
        {
        iter.Set(*e);
        }
      } // next point in region

    } // next element
}

const Element *
Solver::GetElementAtPoint(const VectorType& pt) const
{
  // Add zeros to the end of physical point if necesarry
  Point<Float,MaxGridDimensions> pp;
  for(unsigned int i=0;i<MaxGridDimensions;i++)
    {
    if ( i < pt.size() )
      {
      pp[i]=pt[i];
      }
    else
      {
      pp[i]=0.0;
      }
    }

  InterpolationGridType::IndexType index;

  // Return value only if given point is within the interpolation grid
  if( m_InterpolationGrid->TransformPhysicalPointToIndex(pp,index) )
    {
    //itk::ContinuousIndex<Float,MaxGridDimensions> ci;
    //m_InterpolationGrid->TransformPhysicalPointToContinuousIndex(pp,ci);
    //std::cout<<"In:"<<index<<", ci="<<ci<<", c="<<pp<<"\n";
    return m_InterpolationGrid->GetPixel(index);
    }
  else
    {
    // Return 0, if outside the grid.
    return 0;
    }
}


}} // end namespace itk::fem