itkfemsolvercranknicolson.cxx

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

  unsigned int iter;
  
  Float a,b,d=0.,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;

  Float e=0.0;  // the distance moved on the step before last;

  a=((ax  < cx) ? ax : cx);
  b=((ax  > cx) ? ax : cx);
  
  x=w=v=bx;
  fw=fv=fx=vcl_fabs(EvaluateResidual(x));

  for (iter = 1; iter <=MaxIters; iter++)
  {
    xm=0.5*(a+b);
    tol2=2.0*(tol1=tol*vcl_fabs(x)+ZEPS);
    if (fabs(x-xm) <= (tol2-0.5*(b-a)))
    {
      xmin=x;
      SetEnergyToMin(xmin);
      return fx;
    }

    if (fabs(e) > tol1){
      r=(x-w)*(fx-fv);
      q=(x-v)*(fx-fw);
      p=(x-v)*q-(x-w)*r;
      q=2.0*(q-r);
      if (q>0.0) p = -1.*p;
      q=vcl_fabs(q);
      etemp=e;
      e=d;
      if (fabs(p) >= vcl_fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
         d=CGOLD*(e=(x>=xm ? a-x : b-x));
      else{
        if (q == 0.0) q=q +ZEPS;
        d=p/q;
        u=x+d;
        if (u-a < tol2 || b-u < tol2) d=GSSign(tol1,xm-x); 
      }
  }else{
      d=CGOLD*(e=(x>= xm ? a-x : b-x));
    }
  
    u=(fabs(d) >= tol1 ? x+d : x + GSSign(tol1,d));
    fu=vcl_fabs(EvaluateResidual(u));
    if (fu <= fx){
      if ( u >= x ) a=x; else b=x;
      v=w; w=x;x=u;
      fv=fw; fw=fx; fx=fu;
    } else {
      if (u<x) a = u; else b=u;
      if (fu <= fw || w ==x) {
        v=w;
        w=u;
        fv=fw;
        fw=fu;
      } else if (fu <= fv || v==x || v == w) {
        v=u;
        fv=fu;
      }
    }
  }
  xmin=x;
  SetEnergyToMin(xmin);
  return fx;
}



Element::Float SolverCrankNicolson::GoldenSection(Float tol,unsigned int MaxIters)
{
  // We should now have a, b and c, as well as f(a), f(b), f(c), 
  // where b gives the minimum energy position;

  Float ax, bx, cx;


  FindBracketingTriplet(&ax, &bx, &cx);
  Float xmin,fmin;
  //if (fb!=0.0)
  Float f1,f2,x0,x1,x2,x3;

  Float R=0.6180339;
  Float C=(1.0-R);

  x0=ax;
  x3=cx;
  if (fabs(cx-bx) > vcl_fabs(bx-ax)){
    x1=bx;
    x2=bx+C*(cx-bx);
  } else {
    x2=bx;
    x1=bx-C*(bx-ax);
  }
  f1=vcl_fabs(EvaluateResidual(x1));
  f2=vcl_fabs(EvaluateResidual(x2));
  unsigned int iters=0;
  while (fabs(x3-x0) > tol*(fabs(x1)+vcl_fabs(x2)) && iters < MaxIters)
  {
    iters++;
    if (f2 < f1){
      x0=x1; x1=x2; x2=R*x1+C*x3;
      f1=f2; f2=vcl_fabs(EvaluateResidual(x2));
    } else {
      x3=x2; x2=x1; x1=R*x2+C*x0;
      f2=f1; f1=vcl_fabs(EvaluateResidual(x1));
    }
  }
  if (f1<f2){
    xmin=x1;
    fmin=f1;
  } else {
    xmin=x2;
    fmin=f2;
  }
  
  SetEnergyToMin(xmin);
  return fmin; 
}



void SolverCrankNicolson::SetEnergyToMin(Float xmin)
{
  for (unsigned int j=0; j<NGFN; j++)
    {
    Float SolVal;
    Float FVal;
#ifdef LOCE
    SolVal=xmin*m_ls->GetSolutionValue(j,SolutionTIndex)
      +(1.-xmin)*m_ls->GetSolutionValue(j,SolutionTMinus1Index);   
  
    FVal=xmin*m_ls->GetVectorValue(j,ForceTIndex)
      +(1.-xmin)*m_ls->GetVectorValue(j,ForceTMinus1Index);
#endif
#ifdef TOTE
    SolVal=xmin*m_ls->GetSolutionValue(j,SolutionTIndex);// FOR TOT E
    FVal=xmin*m_ls->GetVectorValue(j,ForceTIndex);
#endif 
    m_ls->SetSolutionValue(j,SolVal,SolutionTIndex);
    m_ls->SetVectorValue(j,FVal,ForceTIndex);
  }

}

Element::Float SolverCrankNicolson::GetDeformationEnergy(Float t)
{
  Float DeformationEnergy=0.0;
  Float iSolVal,jSolVal;

  for (unsigned int i=0; i<NGFN; i++)
    {
// forming  U^T F
#ifdef LOCE
    iSolVal=t*(m_ls->GetSolutionValue(i,SolutionTIndex))
      +(1.-t)*m_ls->GetSolutionValue(i,SolutionTMinus1Index);
#endif
#ifdef TOTE
    iSolVal=t*(m_ls->GetSolutionValue(i,SolutionTIndex))
      +m_ls->GetSolutionValue(i,TotalSolutionIndex);// FOR TOT E
#endif
// forming U^T K U
    Float TempRowVal=0.0;
    for (unsigned int j=0; j<NGFN; j++)
      {
#ifdef LOCE
      jSolVal=t*(m_ls->GetSolutionValue(j,SolutionTIndex))
        +(1.-t)*m_ls->GetSolutionValue(j,SolutionTMinus1Index);
#endif
#ifdef TOTE
      jSolVal=t*(m_ls->GetSolutionValue(j,SolutionTIndex))
        +m_ls->GetSolutionValue(j,TotalSolutionIndex);// FOR TOT E
#endif
      TempRowVal+=m_ls->GetMatrixValue(i,j,SumMatrixIndex)*jSolVal;
      }
    DeformationEnergy+=iSolVal*TempRowVal;
  }
  return DeformationEnergy;
}


Element::Float SolverCrankNicolson::EvaluateResidual(Float t)
{
 
  Float ForceEnergy=0.0,FVal=0.0;
  Float DeformationEnergy=0.0;
  Float iSolVal,jSolVal;

  for (unsigned int i=0; i<NGFN; i++)
  {
// forming  U^T F
#ifdef LOCE
    iSolVal=t*(m_ls->GetSolutionValue(i,SolutionTIndex))
       +(1.-t)*m_ls->GetSolutionValue(i,SolutionTMinus1Index);
    FVal=m_ls->GetVectorValue(i,ForceTIndex);
  FVal=t*FVal+(1.-t)*m_ls->GetVectorValue(i,ForceTMinus1Index);

    ForceEnergy+=iSolVal*FVal;
#endif
#ifdef TOTE
    FVal=FVal+0.0;
    iSolVal=t*(m_ls->GetSolutionValue(i,SolutionTIndex))
        +m_ls->GetSolutionValue(i,TotalSolutionIndex);// FOR TOT E
    ForceEnergy+=iSolVal*(m_ls->GetVectorValue(i,ForceTotalIndex)+
    t*m_ls->GetVectorValue(i,ForceTIndex));// FOR TOT E
#endif
// forming U^T K U
    Float TempRowVal=0.0;
    for (unsigned int j=0; j<NGFN; j++)
    {
#ifdef LOCE
      jSolVal=t*(m_ls->GetSolutionValue(j,SolutionTIndex))
         +(1.-t)*m_ls->GetSolutionValue(j,SolutionTMinus1Index);
#endif
#ifdef TOTE
      jSolVal=t*(m_ls->GetSolutionValue(j,SolutionTIndex))
       +m_ls->GetSolutionValue(j,TotalSolutionIndex);// FOR TOT E
#endif
      TempRowVal+=m_ls->GetMatrixValue(i,j,SumMatrixIndex)*jSolVal;
    }
    DeformationEnergy+=iSolVal*TempRowVal;
  }
  Float Energy=(Float) vcl_fabs(DeformationEnergy-ForceEnergy);
  return Energy;
}

/*
 * Copy solution vector u to the corresponding nodal values, which are
 * stored in node objects). This is standard post processing of the solution.
 */  
void SolverCrankNicolson::AddToDisplacements(Float optimum) 
{
  /*
   * Copy the resulting displacements from 
   * solution vector back to node objects.
   */
  Float maxs=0.0,CurrentTotSolution,CurrentSolution,CurrentForce;
  Float mins2=0.0, maxs2=0.0;
  Float absmax=0.0;

  for(unsigned int i=0;i<NGFN;i++)
  {  
#ifdef TOTE
    CurrentSolution=m_ls->GetSolutionValue(i,SolutionTIndex);
#endif
    if (CurrentSolution < mins2 ){  
      mins2=CurrentSolution;
    }
    else if (CurrentSolution > maxs2 ) {
      maxs2=CurrentSolution;
    } 
    if (fabs(CurrentSolution) > absmax) absmax=vcl_fabs(CurrentSolution);

//  note: set rather than add - i.e. last solution of system not total solution  
#ifdef LOCE
    CurrentSolution=optimum*m_ls->GetSolutionValue(i,SolutionTIndex)
          +(1.-optimum)*m_ls->GetVectorValue(i,SolutionVectorTMinus1Index);
    CurrentForce=optimum*m_ls->GetVectorValue(i,ForceTIndex)
          +(1.-optimum)*m_ls->GetVectorValue(i,ForceTMinus1Index);
    m_ls->SetVectorValue(i,CurrentSolution,SolutionVectorTMinus1Index);   
    m_ls->SetSolutionValue(i,CurrentSolution,SolutionTMinus1Index);   
    m_ls->SetVectorValue(i , CurrentForce, ForceTMinus1Index); // now set t minus one force vector correctly
#endif
#ifdef TOTE
    CurrentSolution=optimum*CurrentSolution;
    CurrentForce=optimum*m_ls->GetVectorValue(i,ForceTIndex);
    m_ls->SetVectorValue(i,CurrentSolution,SolutionVectorTMinus1Index); // FOR TOT E   
    m_ls->SetSolutionValue(i,CurrentSolution,SolutionTMinus1Index);  // FOR TOT E 
    m_ls->SetVectorValue(i,CurrentForce,ForceTMinus1Index);
#endif

    m_ls->AddSolutionValue(i,CurrentSolution,TotalSolutionIndex);
    m_ls->AddVectorValue(i , CurrentForce, ForceTotalIndex);
    CurrentTotSolution=m_ls->GetSolutionValue(i,TotalSolutionIndex);
   
    if ( vcl_fabs(CurrentTotSolution) > maxs ) {
      maxs=vcl_fabs(CurrentTotSolution);
    }

    
  }  
 
  m_CurrentMaxSolution=absmax;
}

/*
 * Compute maximum and minimum solution values.
 */  
void SolverCrankNicolson::PrintMinMaxOfSolution() 
{
  /*
   * Copy the resulting displacements from 
   * solution vector back to node objects.
   */
  Float mins=0.0, maxs=0.0;
  Float mins2=0.0, maxs2=0.0;
  for(unsigned int i=0;i<NGFN;i++)
  {  
    Float CurrentSolution=m_ls->GetSolutionValue(i,SolutionTIndex);
    if (CurrentSolution < mins2 )  mins2=CurrentSolution;
    else if (CurrentSolution > maxs2 )  maxs2=CurrentSolution;
    CurrentSolution=m_ls->GetSolutionValue(i,TotalSolutionIndex);
    if (CurrentSolution < mins )  mins=CurrentSolution;
    else if (CurrentSolution > maxs )  maxs=CurrentSolution;
  }
}

/*
 * Copy solution vector u to the corresponding nodal values, which are
 * stored in node objects). This is standard post processing of the solution.
 */  
void SolverCrankNicolson::AverageLastTwoDisplacements(Float t) 
{
 
  Float maxs=0.0;
  for(unsigned int i=0;i<NGFN;i++)
  {  
    Float temp=m_ls->GetSolutionValue(i,SolutionTIndex);
    Float temp2=m_ls->GetSolutionValue(i,SolutionTMinus1Index);
    Float newsol=t*(temp)+(1.-t)*temp2;
    m_ls->SetSolutionValue(i,newsol,SolutionTMinus1Index);  
    m_ls->SetVectorValue(i,newsol,SolutionVectorTMinus1Index);  
    m_ls->SetSolutionValue(i,newsol,SolutionTIndex);    
    if ( newsol > maxs )  maxs=newsol;
  }  
}

void SolverCrankNicolson::ZeroVector(int which) 
{
  for(unsigned int i=0;i<NGFN;i++)
  {  
    m_ls->SetVectorValue(i,0.0,which);
  }
}

void SolverCrankNicolson::PrintDisplacements() 
{
  std::cout <<  " printing current displacements " << std::endl;
  for(unsigned int i=0;i<NGFN;i++)
  {  
    std::cout << m_ls->GetSolutionValue(i,TotalSolutionIndex) << std::endl;
  }
}

void SolverCrankNicolson::PrintForce() 
{
  std::cout <<  " printing current forces " << std::endl;
  for(unsigned int i=0;i<NGFN;i++)
  {  
    std::cout << m_ls->GetVectorValue(i,ForceTIndex) << std::endl;
  }
}




}} // end namespace itk::fem