Cmesh_generate_regularFem.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 CMESH_GENERATE_REGULARFEM
#define CMESH_GENERATE_REGULARFEM

#include "Cmesh_generate.h"

template<typename T, typename Timg> 
class Cmesh_generate_regularFem : public Cmesh_generate<T,Timg>
{

  
  public:
 
  Cmesh_generate_regularFem() : Cmesh_generate<T,Timg>()
  {// constructor
    
    this->class_name = "Cmesh_generate_regularFem";
    
  }
  
  virtual int get_attributes(CParameterNetCDF &fp)
  {// requires to load pitch[x,y,z]
  
    // load components
    Cmesh<T,Timg>::get_attributes(fp);
    
    int error(0);
    std::string att;

    // store pitch sizes within vector[3] named m_pitch. The last value is one in the 2D case.
    att = "pitch";
    error = fp.loadAttribute(att,this->m_pitch);
    for (int i=this->m_pitch.size(); i<3; i++){this->m_pitch.push_back(1);}
    this->m_win.assign(3,0);
    for (int i=0; i<3; i++){this->m_win[i] = this->m_pitch[i];}
    
    if (this->m_verbose)
    {printf("\tMesh resolution : [%i %i %i]\n",this->m_pitch[0], this->m_pitch[1], this->m_pitch[2]);}

  }
  
  virtual int assign(const Cimage<Timg> &oImage)
  {// initialize attribute dimensions

    //grid dimensions
    this->m_grid.assign(3,1);
    for (int i=0; i<this->num_var(); i++)
    {
      this->m_grid[i] = (int)ceil((float)oImage.cur_dim(i)/(float)this->m_pitch[i])+1;
    }

    //degree of freedom
    this->m_dof = this->m_grid[0]*this->m_grid[1]*this->m_grid[2]*this->num_var();
    if (this->m_verbose){printf("\tDegree Of Freedom : %i\n",this->m_dof);}
    
    //node table
    this->m_node.assign(this->num_var(),this->m_grid[0],this->m_grid[1],this->m_grid[2],1,0);
    if (this->m_verbose){printf("\tNode table : [%i](%i,%i,%i)\n",this->m_node.size(), this->m_node[0].width(), this->m_node[0].height(), this->m_node[0].depth());}
    
    //connectivity table
    if (this->num_var() == 2)
    {//2D
      this->m_connect.assign((this->m_grid[0]-1)*(this->m_grid[1]-1),2,2,1,1,0);
    }
    else
    {//3D
      this->m_connect.assign((this->m_grid[0]-1)*(this->m_grid[1]-1)*(this->m_grid[2]-1),2,2,2,1,0);
    }
    if (this->m_verbose){printf("\tElement table : [%i](%i,%i,%i)\n",this->m_connect.size(), this->m_connect[0].width(), this->m_connect[0].height(), this->m_connect[0].depth());}
    
    //number of pixels
    this->m_numPix = oImage.m_curImg[0].size();
    
    return 0;

  }
  
  virtual void set_nodes(const Cimage<Timg> &oImage, CImgList<T> &nodes)
  {// node table definition
  
    std::vector<int> x;
    x.assign(3,0);
    
    #pragma omp for schedule(runtime)
    for (int nd=0; nd<nodes[0].size(); nd++)
    {//fill node table
      
      //find back spatial coordonnates of nodes
      nodes[0].contains(nodes[0][nd],x[0],x[1],x[2]);

      cimglist_for(nodes, dim)
      {//allow different element size at image borders to obtain nodes at each image corner
      
        if (x[dim] == this->m_grid[dim]-1)
        {nodes[dim][nd] = oImage.cur_dim(dim)-1;}
        else
        {nodes[dim][nd] = x[dim] * this->m_pitch[dim];}
        
      }

    }
    

  }
  
  virtual void set_connectivity(const Cimage<Timg> &oImage, const CImgList<T> &nodes, CImgList<int> &elements)
  {// element and connectivity table definition
    
    std::vector<int> x;
    x.assign(3,0);
    
    #pragma omp for schedule(runtime)
    for (int elem=0; elem<elements.size(); elem++)
    {//fill element table
      
      //find back spatial coordonnates within element table
      x[0] = elem%(this->m_grid[0]-1);
      x[1] = (elem/(this->m_grid[0]-1))%(this->m_grid[1]-1);
      x[2] = elem/((this->m_grid[1]-1)*(this->m_grid[0]-1));
      
      //fond associated nodes
      cimg_forZ(elements[elem], c)
      {
        elements(elem, 0,0,c) = nodes[0].offset(x[0],x[1],x[2]+c,0);
        elements(elem, 1,0,c) = nodes[0].offset(x[0]+1,x[1],x[2]+c,0);
        elements(elem, 1,1,c) = nodes[0].offset(x[0]+1,x[1]+1,x[2]+c,0);
        elements(elem, 0,1,c) = nodes[0].offset(x[0],x[1]+1,x[2]+c,0);
      }
      
    }
    
    

  }

  virtual int cropList(const Cimage<Timg> &oImage, const int &elem, CImg<Timg> &list, std::vector<T> &coord0) const
  {//extract sub_box from global image one using node and element table
  
    std::vector<T> coord1;
    coord0.assign(3,0);coord1.assign(3,0);
    int msk(1);

    cimglist_for(this->m_node, dim)
    {//set corner box positions
      
      //upper left and lower righ nodes
      coord0[dim] = this->m_node[dim][this->m_connect[elem].front()];
      coord1[dim] = this->m_node[dim][this->m_connect[elem].back()]-1;
      
      if (coord0[dim] == 0)
      {coord0[dim]++;}
      
    }
    
    double unmasked((coord1[0]-coord0[0]+1)*(coord1[1]-coord0[1]+1)*(coord1[2]-coord0[2]+1));
    
    if (!oImage.m_curMask.is_empty())
    {//number of unmasked pixels within element
      unmasked = oImage.m_curMask.get_crop(coord0[0],coord0[1],coord0[2],0,coord1[0],coord1[1],coord1[2],0,0).sum();
    }
    
    if (unmasked >= 10*(coord1[0]-coord0[0]+1)*(coord1[1]-coord0[1]+1)*(coord1[2]-coord0[2]+1)/100) //(coord1[0]-coord0[0]+1)*(coord1[1]-coord0[1]+1)*(coord1[2]-coord0[2]+1))
    {list.assign(coord1[0]-coord0[0]+1, coord1[1]-coord0[1]+1, coord1[2]-coord0[2]+1, 1, 0);}
    else
    {msk = 0;}

    return msk;
    
  }

  virtual int grid_dims(const int &i) const
  {
  
    switch(i)
    {
      case 0: return this->m_grid[0]; break;
      case 1: return this->m_grid[1]; break;
      case 2: return this->m_grid[2]; break;
      case 3: return 1; break;
      default:
      {
        std::cerr<<this->class_name<<"::"<<__func__<<": error: id="<<i<<" which is over [0-3]\n"<<std::flush;
        return 0;
      }
    }
  
  }
 
};

#endif