amroc/amr/GFMExactSolution.h

Go to the documentation of this file.

// -*- C++ -*-

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

#ifndef AMROC_GFM_EXACT_SOLUTION_H
#define AMROC_GFM_EXACT_SOLUTION_H

#include "ExactSolution.h"

template <class VectorType, class FixupType, class FlagType, int dim>
class GFMExactSolution : public ExactSolution<VectorType,dim> {
  typedef ExactSolution<VectorType,dim> base;
  typedef typename VectorType::InternalDataType DataType;

public:
  typedef typename base::vec_grid_fct_type vec_grid_fct_type;  
  typedef typename base::grid_fct_type grid_fct_type;
  typedef typename base::grid_data_type grid_data_type;

  typedef AMRGFMSolver<VectorType,FixupType,FlagType,dim> gfm_solver_type;
  typedef GhostFluidMethod<VectorType,dim> gfm_type;

  GFMExactSolution(gfm_solver_type& solver) : base(), _solver(solver) {} 

  virtual ~GFMExactSolution() {} 

  virtual void CalculateNorm(vec_grid_fct_type& u, double*& Error) {
    if (NGFM()>0) 
      for (int Level=0; Level<=FineLevel(base::GH()); Level++) {
        int Time = CurrentTime(base::GH(),Level); 
        int TStep = TimeStep(u,Level);      
        forall (u,Time+TStep,Level,c) 
          for (int g=0; g<NGFM(); g++) {
            if (!GFM(g).IsUsed()) continue;
            if (base::Dim() == 1) {
              BeginFastIndex1(phi, Phi(g)(Time,Level,c).bbox(), 
                            Phi(g)(Time,Level,c).data(), DataType);
              BeginFastIndex1(u, u(Time+TStep,Level,c).bbox(), 
                              u(Time+TStep,Level,c).data(), VectorType);
              for_1 (n, Phi(g)(Time,Level,c).bbox(), Phi(g)(Time,Level,c).bbox().stepsize())
                if (FastIndex1(phi,n)<-Cutoff(g)) 
                  FastIndex1(u,n) = static_cast<VectorType>(0.0);
              end_for
              EndFastIndex1(u);
              EndFastIndex1(phi);
            }      
            else if (base::Dim() == 2) {
              BeginFastIndex2(phi, Phi(g)(Time,Level,c).bbox(), 
                              Phi(g)(Time,Level,c).data(), DataType);
              BeginFastIndex2(u, u(Time+TStep,Level,c).bbox(), 
                              u(Time+TStep,Level,c).data(), VectorType);
              for_2 (n, m, Phi(g)(Time,Level,c).bbox(), Phi(g)(Time,Level,c).bbox().stepsize())
                if (FastIndex2(phi,n,m)<-Cutoff(g)) 
                  FastIndex2(u,n,m) = static_cast<VectorType>(0.0);
              end_for
              EndFastIndex2(u);
              EndFastIndex2(phi);
            }
            else if (base::Dim() == 3) {
              BeginFastIndex3(phi, Phi(g)(Time,Level,c).bbox(), 
                              Phi(g)(Time,Level,c).data(), DataType);
              BeginFastIndex3(u, u(Time+TStep,Level,c).bbox(), 
                              u(Time+TStep,Level,c).data(), VectorType);
              for_3 (n, m, l, Phi(g)(Time,Level,c).bbox(), Phi(g)(Time,Level,c).bbox().stepsize())
                if (FastIndex3(phi,n,m,l)<-Cutoff(g)) 
                  FastIndex3(u,n,m,l) = static_cast<VectorType>(0.0);
              end_for
              EndFastIndex3(u);
              EndFastIndex3(phi);
            }
          }
        end_forall
      }
    base::CalculateNorm(u,Error);
  }

  virtual void CalculateNorm(grid_fct_type& u, double*& Error) {
    if (NGFM()>0) 
      for (int Level=0; Level<=FineLevel(base::GH()); Level++) {
        int Time = CurrentTime(base::GH(),Level); 
        int TStep = TimeStep(u,Level);      
        forall (u,Time+TStep,Level,c) 
          for (int g=0; g<NGFM(); g++) {
            if (!GFM(g).IsUsed()) continue;
            if (base::Dim() == 1) {
              BeginFastIndex1(phi, Phi(g)(Time,Level,c).bbox(), 
                            Phi(g)(Time,Level,c).data(), DataType);
              BeginFastIndex1(u, u(Time+TStep,Level,c).bbox(), 
                              u(Time+TStep,Level,c).data(), DataType);
              for_1 (n, Phi(g)(Time,Level,c).bbox(), Phi(g)(Time,Level,c).bbox().stepsize())
                if (FastIndex1(phi,n)<-Cutoff(g)) 
                  FastIndex1(u,n) = static_cast<DataType>(0.0);
              end_for
              EndFastIndex1(u);
              EndFastIndex1(phi);
            }      
            else if (base::Dim() == 2) {
              BeginFastIndex2(phi, Phi(g)(Time,Level,c).bbox(), 
                              Phi(g)(Time,Level,c).data(), DataType);
              BeginFastIndex2(u, u(Time+TStep,Level,c).bbox(), 
                              u(Time+TStep,Level,c).data(), DataType);
              for_2 (n, m, Phi(g)(Time,Level,c).bbox(), Phi(g)(Time,Level,c).bbox().stepsize())
                if (FastIndex2(phi,n,m)<-Cutoff(g)) 
                  FastIndex2(u,n,m) = static_cast<DataType>(0.0);
              end_for
              EndFastIndex2(u);
              EndFastIndex2(phi);
            }
            else if (base::Dim() == 3) {
              BeginFastIndex3(phi, Phi(g)(Time,Level,c).bbox(), 
                              Phi(g)(Time,Level,c).data(), DataType);
              BeginFastIndex3(u, u(Time+TStep,Level,c).bbox(), 
                              u(Time+TStep,Level,c).data(), DataType);
              for_3 (n, m, l, Phi(g)(Time,Level,c).bbox(), Phi(g)(Time,Level,c).bbox().stepsize())
                if (FastIndex3(phi,n,m,l)<-Cutoff(g)) 
                  FastIndex3(u,n,m,l) = static_cast<DataType>(0.0);
              end_for
              EndFastIndex3(u);
              EndFastIndex3(phi);
            }
          }
        end_forall
      }
    base::CalculateNorm(u,Error);
  }

  inline gfm_type& GFM(const int& n) { return _solver.GFM(n); }
  inline const int& NGFM() { return _solver.NGFM(); }
  inline grid_fct_type& Phi(const int& n) { return GFM(n).Phi(); }
  inline const DataType& Cutoff(const int& n) { return GFM(n).Boundary().Cutoff(); }
protected:
  gfm_solver_type& _solver;
};

#endif

---