Ccorrelation_opticalFlow_fem_newton.h

Go to the documentation of this file.

/**********************************************************************
 * Copyright (C) 2012, The YaDICs Project Developers. 
 * See the COPYRIGHT file at the top-level directory of this distribution ./COPYRIGHT.
 * See ./COPYING file for copying and redistribution conditions.
 * 
 * This file is part of YaDICs.
 * YaDICs is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * YaDICs is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with YaDICs.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Information about how to use the software are provided at http://yadic.univ-lille1.fr/
 **********************************************************************/


#ifndef CCORRELATION_OPTICALFLOW_FEM_NEWTON
#define CCORRELATION_OPTICALFLOW_FEM_NEWTON


#include "Ccorrelation_opticalFlow_fem.h"



template<typename T,typename Timg>
class Ccorrelation_opticalFlow_fem_newton: public Ccorrelation_opticalFlow_fem<T,Timg>
{
  
  public:

    Ccorrelation_opticalFlow_fem_newton() : Ccorrelation_opticalFlow_fem<T,Timg>()
    {// constructor
      
      this->class_name = "Correlation : Optical Flow Finite Element Method using Newton strategy";

    }
      
    virtual void solve(CImgList<T> &K, const Cmesh<T,Timg> &oMesh, CImgList<T> &solution)
    {// K.U = F

      std::vector<int> count;
      CImgList<T> Kc(2);
      count.clear();
      
      for (int nd = 0; nd<oMesh.m_node[0].size(); nd++)
      {
        if (K[1][oMesh.num_var()*nd] != 0)
        {
          count.push_back(nd);
        }
      }

      if (count.size()*oMesh.num_var() != K[0].width())
      {
        if(this->m_verbose)
        {
          printf("Stiffness tensor densification required (because of mask)\n");
          std::cout<<"Densification : from "<<count.size()*oMesh.num_var()<<" to "<<K[0].width()<<" \n";
        }
      
        //assign compact k and f matrix
        Kc[0].assign(count.size()*oMesh.num_var(), count.size()*oMesh.num_var(), 1, 1, 0);
        Kc[1].assign(1, count.size()*oMesh.num_var(), 1, 1, 0);
      
        for (int dim=0; dim<oMesh.num_var(); dim++)
        {
          for (int i=0; i<count.size(); i++)
          {//transfert porous K and F matrix data within compact k and f ones
            for (int j=0; j<count.size(); j++)
            {
              Kc[0](i*oMesh.num_var()+dim,j*oMesh.num_var()+dim) = K[0](count[i]*oMesh.num_var()+dim, count[j]*oMesh.num_var()+dim);
            }
            Kc[1][i*oMesh.num_var()+dim] = K[1][count[i]*oMesh.num_var()+dim];
          }
        }
        
      }
      else
      {
        Kc = K.get_shared();
      }

      Kc[1].solve(Kc[0]);
       
      for (int node=0; node<count.size(); node++)
      {//put solution vector in oField.m_mode[dim][node] framework
        for (int dim=0; dim<oMesh.num_var(); dim++)
        {
          solution[dim][count[node]] += Kc[1][node*oMesh.num_var()+dim];
        }
      }
     
    }  

    virtual void assignG(const Cmesh<T,Timg> &oMesh, CImgList<T> &K, CImgList<T> &Ke)
    {//assign approriate size for Ke and Fe tensors

        Ccorrelation_opticalFlow<T,Timg>::assignG(oMesh, K, Ke);
        
        K[0].assign(oMesh.num_node()*oMesh.num_var(),oMesh.num_node()*oMesh.num_var(),1,1, 0);
        K[1].assign(1,oMesh.num_node()*oMesh.num_var(),1,1, 0);
        
        Ke[0].assign(oMesh.num_mod(),oMesh.num_mod(),oMesh.num_var(),this->m_threadNB, 0);
        Ke[1].assign(1,oMesh.num_mod(),oMesh.num_var(),this->m_threadNB, 0);
   
    }

    virtual void assembly(const int &thread, const Cmesh<T,Timg> &oMesh, const int &elem, const CImgList<T> &Ke, CImgList<T> &K)
    {//perform the assembly of K and F

      //Ke[0] = Ke[modi, modj, dims]
      //Ke[1] = Fe[1, modi, dims]
    
      cimg_forZ(Ke[0], dim)
      {
        cimg_forX(Ke[0], modi)
        {
          cimg_forY(Ke[0], modj)
          {
            #pragma omp critical
            {
              K[0](oMesh.m_connect[elem][modi]*Ke[0].depth()+dim,oMesh.m_connect[elem][modj]*Ke[0].depth()+dim) += Ke[0](modi,modj,dim,thread);
            }
          }

          #pragma omp critical
          {
            K[1][oMesh.m_connect[elem][modi]*Ke[0].depth()+dim] += Ke[1](0,modi,dim,thread);
          }
        }
      }
      
    }
  
  
};
 
#endif