amroc/converters/Evaluators/NormEvaluator.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2002 Ralf Deiterding
// Brandenburgische Universitaet Cottbus
//
// Copyright (C) 2003-2007 California Institute of Technology
// Ralf Deiterding, ralf@cacr.caltech.edu

// Copyright (C) 2010 Oak Ridge National Laboratory
// Ralf Deiterding, deiterdingr@ornl.gov

#ifndef AMROC_NORMEVALUATOR_H
#define AMROC_NORMEVALUATOR_H

#include <iostream>
#include <cstdio>
#include <cfloat>

#include "Evaluator.h"

template <class VectorType>
class NormEvaluator : public Evaluator<VectorType> {
  typedef Evaluator<VectorType> base;
  typedef typename VectorType::InternalDataType DataType;

public:
  typedef typename base::vector_block_type vector_block_type;
  typedef typename base::data_block_type data_block_type;

  NormEvaluator() : base() {
    std::sprintf(base::title,"Norm evaluation of scalar functions");
  }

  virtual void SetUp(VisualGridBase *gr) {
    EvaluatorBase::SetUp(gr);

    std::sprintf(&(base::tkeys[0*LEN_TKEYS]), "Distribution            ");
    std::sprintf(&(base::tkeys[1*LEN_TKEYS]), "Levels                  ");
    std::sprintf(&(base::tkeys[2*LEN_TKEYS]), "|F1-F2|                 ");
    std::sprintf(&(base::tkeys[3*LEN_TKEYS]), "|F1-F2| dx dy dz        ");
    std::sprintf(&(base::tkeys[4*LEN_TKEYS]), "|F1-F2|^2 dx dy dz      ");
    base::fkeys[0] = base::Grid->FKeys();
    base::fkeys[1] = base::Grid->FKeys();
    base::fkeys[2] = base::Grid->FKeys();
    base::fkeys[3] = base::Grid->FKeys();
    base::fkeys[4] = base::Grid->FKeys();
    base::ikeys[0] = 105; base::ikeys[1] = 115; 
    base::ikeys[2] = 109; base::ikeys[3] = 49; 
    base::ikeys[4] = 50; 
  }    

  virtual int NKeys() const { return 5; }

  virtual void SetScalCells(int key,float v[],int& offset, 
                            vector_block_type& DB, float* dx, bool* CompRead) {
    int idx=0;
    int keym1 = key-1;
    if (key>=3 && key<=5) {
      BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
      float f1, f2;
      float fac = dx[0]*DB.stepsize(0)*dx[1]*DB.stepsize(1)*dx[2]*DB.stepsize(2);
      for_3 (i,j,k,DB.bbox(),DB.bbox().stepsize())
        v[idx+offset]= FastIndex3(dat,i,j,k)(0);
        f1 = FastIndex3(dat,i,j,k)(0);
        f2 = FastIndex3(dat,i,j,k)(1);
        if (f1<FLT_MAX && f2<FLT_MAX) {
          if (key==3)
            v[idx+offset] = std::fabs(f1-f2);
          else if (key==4)
            v[idx+offset] = std::fabs(f1-f2)*fac;
          else if (key==5)
            v[idx+offset] = (f1-f2)*(f1-f2)*fac;
        }
        base::mvals[2*keym1]   = Min(base::mvals[2*keym1],   v[idx+offset]);
        base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
        idx++;        
      end_for
      EndFastIndex3(dat); 
    }
    offset+=idx;
  }
 
  virtual void SetScalNodes(int key,float v[],int& offset,
                            vector_block_type& DB, const BBox& bbox, 
                            float* dx, bool* CompRead) {
    int idx=0;
    int keym1 = key-1;
    float fac = dx[0]*DB.stepsize(0)*dx[1]*DB.stepsize(1)*dx[2]*DB.stepsize(2);
    BeginFastIndex3(dat, DB.bbox(), DB.data(), VectorType);
    const Coords& step = bbox.stepsize();
    switch(key) {
    case 3: 
      {
        for_3 (i,j,k,bbox,step)
          BBox AddOver(3,i-step(0),j-step(1),k-step(2),
                       i,j,k,step(0),step(1),step(2));
          AddOver *= DB.bbox();
          
          v[idx+offset] = 0.0;
          for_3 (l, m, n, AddOver, step)
            if (FastIndex3(dat,l,m,n)(0)<FLT_MAX && FastIndex3(dat,l,m,n)(1)<FLT_MAX) {
              v[idx+offset] = Max(std::fabs(FastIndex3(dat,l,m,n)(0)-FastIndex3(dat,l,m,n)(1)),v[idx+offset]);
            }
          end_for

          base::mvals[2*keym1]   = Min(base::mvals[2*keym1],   v[idx+offset]);
          base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
          idx++;              
        end_for 
        break;
      }

    case 4:
      {
        for_3 (i,j,k,bbox,step)
          BBox AddOver(3,i-step(0),j-step(1),k-step(2),
                       i,j,k,step(0),step(1),step(2));
          AddOver *= DB.bbox();
        
          float c=0;
          v[idx+offset] = 0.0;
          for_3 (l, m, n, AddOver, step)
            if (FastIndex3(dat,l,m,n)(0)<FLT_MAX && FastIndex3(dat,l,m,n)(1)<FLT_MAX) {
              v[idx+offset] += std::fabs(FastIndex3(dat,l,m,n)(0)-FastIndex3(dat,l,m,n)(1));
              c += 1.0;
            }
          end_for
          if (c>0) {
            v[idx+offset] /= c; v[idx+offset] *= fac;
          }
          else v[idx+offset] = FLT_MAX;
        
          base::mvals[2*keym1]   = Min(base::mvals[2*keym1],   v[idx+offset]);
          base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
          idx++;              
        end_for 
        break;
      }

    case 5:
      {
        for_3 (i,j,k,bbox,step)
          BBox AddOver(3,i-step(0),j-step(1),k-step(2),
                       i,j,k,step(0),step(1),step(2));
          AddOver *= DB.bbox();

          float c=0;
          v[idx+offset] = 0.0;
          for_3 (l, m, n, AddOver, step)
            if (FastIndex3(dat,l,m,n)(0)<FLT_MAX && FastIndex3(dat,l,m,n)(1)<FLT_MAX) {
              v[idx+offset] += std::fabs(FastIndex3(dat,l,m,n)(0)-FastIndex3(dat,l,m,n)(1));
              v[idx+offset] *= v[idx+offset];
              c += 1.0;
            }
          end_for
          if (c>0) {
            v[idx+offset] /= c; v[idx+offset] *= fac;
          }
          else v[idx+offset] = FLT_MAX;
        
          base::mvals[2*keym1]   = Min(base::mvals[2*keym1],   v[idx+offset]);
          base::mvals[2*keym1+1] = Max(base::mvals[2*keym1+1], v[idx+offset]);
          idx++;              
        end_for 
        break;
      }
    default:
      break;
    }
    EndFastIndex3(dat); 
    offset+=idx;
  }
};


#endif

---