ScdElementData.hpp

Go to the documentation of this file.

#ifndef SCD_ELEMENT_DATA_HPP
#define SCD_ELEMENT_DATA_HPP

//
// Class: ScdElementData
//
// Purpose: represent a rectangular element of mesh
//
// A ScdElement represents a rectangular element of mesh, including both vertices and
// elements, and the parametric space used to address that element.  Vertex data,
// i.e. coordinates, may not be stored directly in the element, but the element returns
// information about the vertex handles of vertices in the element.  Vertex and element
// handles associated with the element are each contiguous.

#include "SequenceData.hpp"
#include "moab/HomXform.hpp"
#include "moab/CN.hpp"
#include "ScdVertexData.hpp"
#include "Internals.hpp"
#include "moab/Range.hpp"

#include <vector>
#include <algorithm>

namespace moab {

class ScdElementData : public SequenceData
{

class VertexDataRef
{
private:
  HomCoord minmax[2];
  HomXform xform, invXform;
  ScdVertexData *srcSeq;
public:
  friend class ScdElementData;
  
  VertexDataRef(const HomCoord &min, const HomCoord &max,
                const HomXform &tmp_xform, ScdVertexData *this_seq);
    
  bool contains(const HomCoord &coords) const;
};

private:

  HomCoord boxParams[3];

  int dIJK[3];
  
  int dIJKm1[3];

  int isPeriodic[2];
  
  ScdElementData();

  std::vector<VertexDataRef> vertexSeqRefs;

public:

  ScdElementData( EntityHandle start_handle,
                  const int imin, const int jmin, const int kmin,
                  const int imax, const int jmax, const int kmax,
                  int *is_periodic);
  
  virtual ~ScdElementData();
  
  inline EntityHandle get_vertex(const HomCoord &coords) const;
  inline EntityHandle get_vertex(int i, int j, int k) const
    { return get_vertex(HomCoord(i,j,k)); }
  
  EntityHandle get_element(const int i, const int j, const int k) const;
  
  const HomCoord &min_params() const;

  const HomCoord &max_params() const;
  
  void param_extents(int &di, int &dj, int &dk) const;

  ErrorCode get_params(const EntityHandle ehandle,
                          int &i, int &j, int &k) const;

  inline int is_periodic_i() const {return isPeriodic[0];};
  
  inline int is_periodic_j() const {return isPeriodic[1];};

  inline void is_periodic(int is_p[2]) const 
      {is_p[0] = isPeriodic[0]; is_p[1] = isPeriodic[1];}
  
  int i_min() const {return (boxParams[0].hom_coord())[0];}
  int j_min() const {return (boxParams[0].hom_coord())[1];}
  int k_min() const {return (boxParams[0].hom_coord())[2];}
  int i_max() const {return (boxParams[1].hom_coord())[0];}
  int j_max() const {return (boxParams[1].hom_coord())[1];}
  int k_max() const {return (boxParams[1].hom_coord())[2];}

  bool boundary_complete() const;

  inline bool contains(const HomCoord &coords) const;

  inline bool contains_vertex(const HomCoord &coords) const;

  inline ErrorCode get_params_connectivity(const int i, const int j, const int k,
                                       std::vector<EntityHandle>& connectivity) const;
  
  ErrorCode add_vsequence(ScdVertexData *vseq, 
                             const HomCoord &p1, const HomCoord &q1,
                             const HomCoord &p2, const HomCoord &q2,
                             const HomCoord &p3, const HomCoord &q3,
                             bool bb_input = false,
                             const HomCoord &bb_min = HomCoord::unitv[0],
                             const HomCoord &bb_max = HomCoord::unitv[0]);


  SequenceData* subset( EntityHandle start, 
                        EntityHandle end,
                        const int* sequence_data_sizes,
                        const int* tag_data_sizes ) const;
  
  static EntityID calc_num_entities( EntityHandle start_handle,
                                       int irange,
                                       int jrange,
                                     int krange,
                                     int *is_periodic = NULL);

  unsigned long get_memory_use() const;
};

inline EntityHandle ScdElementData::get_element(const int i, const int j, const int k) const
{
  return start_handle() + (i-i_min()) + (j-j_min())*dIJKm1[0] + (k-k_min())*dIJKm1[0]*dIJKm1[1];
}

inline const HomCoord &ScdElementData::min_params() const
{
  return boxParams[0];
}

inline const HomCoord &ScdElementData::max_params() const
{
  return boxParams[1];
}

inline void ScdElementData::param_extents(int &di, int &dj, int &dk) const
{
  di = dIJK[0];
  dj = dIJK[1];
  dk = dIJK[2];
}

inline ErrorCode ScdElementData::get_params(const EntityHandle ehandle,
                                              int &i, int &j, int &k) const
{
  if (TYPE_FROM_HANDLE(ehandle) != TYPE_FROM_HANDLE(start_handle())) return MB_FAILURE;

  int hdiff = ehandle - start_handle();

    // use double ?: test below because on some platforms, both sides of the : are
    // evaluated, and if dIJKm1[1] is zero, that'll generate a divide-by-zero
  k = (dIJKm1[1] > 0 ? hdiff / (dIJKm1[1] > 0 ? dIJKm1[0]*dIJKm1[1] : 1) : 0);
  j = (hdiff - (k*dIJKm1[0]*dIJKm1[1])) / dIJKm1[0];
  i = hdiff % dIJKm1[0];

  k += boxParams[0].k();
  j += boxParams[0].j();
  i += boxParams[0].i();

  return (ehandle >= start_handle() &&
          ehandle < start_handle()+size() &&
          i >= i_min() && i <= i_max() &&
          j >= j_min() && j <= j_max() &&
          k >= k_min() && k <= k_max()) ? MB_SUCCESS : MB_FAILURE;
}

inline bool ScdElementData::contains(const HomCoord &temp) const 
{
    // upper bound is < instead of <= because element params max is one less
    // than vertex params max, except in case of 2d or 1d sequence
  return ((dIJKm1[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i()+dIJKm1[0]) &&
          ((!dIJKm1[1] && temp.j() == boxParams[1].j()) || 
           (dIJKm1[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j()+dIJKm1[1])) &&
          ((!dIJKm1[2] && temp.k() == boxParams[1].k()) || 
           (dIJKm1[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k()+dIJKm1[2])));
}
  
inline bool ScdElementData::contains_vertex(const HomCoord &temp) const 
{
    // upper bound is < instead of <= because element params max is one less
    // than vertex params max, except in case of 2d or 1d sequence
  return ((dIJK[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i()+dIJK[0]) &&
          ((!dIJK[1] && temp.j() == boxParams[1].j()) || 
           (dIJK[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j()+dIJK[1])) &&
          ((!dIJK[2] && temp.k() == boxParams[1].k()) || 
           (dIJK[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k()+dIJK[2])));
}
  
inline bool ScdElementData::VertexDataRef::contains(const HomCoord &coords) const 
{
  return (minmax[0] <= coords && minmax[1] >= coords);
}

inline ScdElementData::VertexDataRef::VertexDataRef(const HomCoord &this_min, const HomCoord &this_max,
                                    const HomXform &tmp_xform, ScdVertexData *this_seq)
    : xform(tmp_xform), invXform(tmp_xform.inverse()), srcSeq(this_seq)
{
  minmax[0] = HomCoord(this_min);
  minmax[1] = HomCoord(this_max); 
}

inline EntityHandle ScdElementData::get_vertex(const HomCoord &coords) const
{
  assert(boundary_complete());
   for (std::vector<VertexDataRef>::const_iterator it = vertexSeqRefs.begin();
        it != vertexSeqRefs.end(); it++) {
     if ((*it).minmax[0] <= coords && (*it).minmax[1] >= coords) {
         // first get the vertex block-local parameters
       HomCoord local_coords = coords / (*it).xform;

       assert((*it).srcSeq->contains(local_coords));

      // now get the vertex handle for those coords
       return (*it).srcSeq->get_vertex(local_coords);
     }
   }
   
     // got here, it's an error
   return 0;
}

inline ErrorCode ScdElementData::add_vsequence(ScdVertexData *vseq, 
                                                 const HomCoord &p1, const HomCoord &q1,
                                                 const HomCoord &p2, const HomCoord &q2, 
                                                 const HomCoord &p3, const HomCoord &q3,
                                                 bool bb_input,
                                                 const HomCoord &bb_min,
                                                 const HomCoord &bb_max)
{
    // compute the transform given the vseq-local parameters and the mapping to
    // this element sequence's parameters passed in minmax
  HomXform M;
  M.three_pt_xform(p1, q1, p2, q2, p3, q3);
  
    // min and max in element seq's parameter system may not be same as those in 
    // vseq's system, so need to take min/max

  HomCoord minmax[2];
  if (bb_input) {
    minmax[0] = bb_min;
    minmax[1] = bb_max;
  }
  else {
    minmax[0] = vseq->min_params() * M;
    minmax[1] = vseq->max_params() * M;
  }
  
    // check against other vseq's to make sure they don't overlap
  for (std::vector<VertexDataRef>::const_iterator vsit = vertexSeqRefs.begin();
       vsit != vertexSeqRefs.end(); vsit++) 
    if ((*vsit).contains(minmax[0]) || (*vsit).contains(minmax[1])) 
      return MB_FAILURE;
    
  HomCoord tmp_min(std::min(minmax[0].i(), minmax[1].i()), 
                   std::min(minmax[0].j(), minmax[1].j()), 
                   std::min(minmax[0].k(), minmax[1].k()));
  HomCoord tmp_max(std::max(minmax[0].i(), minmax[1].i()), 
                   std::max(minmax[0].j(), minmax[1].j()), 
                   std::max(minmax[0].k(), minmax[1].k()));

  
    // set up a new vertex sequence reference
  VertexDataRef tmp_seq_ref(tmp_min, tmp_max, M, vseq);

    // add to the list
  vertexSeqRefs.push_back(tmp_seq_ref);
  
  return MB_SUCCESS;
}

inline ErrorCode ScdElementData::get_params_connectivity(const int i, const int j, const int k,
                                                           std::vector<EntityHandle>& connectivity) const
{
  if (contains(HomCoord(i, j, k)) == false) return MB_FAILURE;
  
  int ip1 = (isPeriodic[0] ? (i+1)%dIJKm1[0] : i+1),
      jp1 = (isPeriodic[1] ? (j+1)%dIJKm1[1] : j+1);
  
  connectivity.push_back(get_vertex(i, j, k));
  connectivity.push_back(get_vertex(ip1, j, k));
  if (CN::Dimension(TYPE_FROM_HANDLE(start_handle())) < 2) return MB_SUCCESS;
  connectivity.push_back(get_vertex(ip1, jp1, k));
  connectivity.push_back(get_vertex(i, jp1, k));
  if (CN::Dimension(TYPE_FROM_HANDLE(start_handle())) < 3) return MB_SUCCESS;
  connectivity.push_back(get_vertex(i, j, k+1));
  connectivity.push_back(get_vertex(ip1, j, k+1));
  connectivity.push_back(get_vertex(ip1, jp1, k+1));
  connectivity.push_back(get_vertex(i, jp1, k+1));
  return MB_SUCCESS;
}
  
} // namespace moab

#endif