weno/applications/euler/2d/Angle/src/Problem.h

// -*- C++ -*-
// vim: syntax=cpp

// $Id: Problem.h,v 1.4 2007/04/25 20:42:47 ferrante Exp $

// To make sure that the tip of the wedge remains smooth when the grid size is
// changed, the variable TipRadius is defined as
//
//                 radius of the tip
//    TipRadius = -------------------
//                  MIN ( Delta X )
//
// aferrante 28mar07 - found that as implemented by mhesseli TipRadius is
//                     TipRadius = radius of the tip * MIN ( Delta X)
// but I want TipRadius to be a parameter specified in input independent on the mesh
// thus I remove (*h) in the formula of TipRadius below

#ifndef AMROC_PROBLEM_H
#define AMROC_PROBLEM_H

#define DIM                            2
#define NPATCHES                       1

#include <cfloat>
#include <sstream>

#include "WENOProblem.h"

#define OWN_GFMAMRSOLVER
#define OWN_FLAGGING

#include "WENOStdGFMProblem.h"

#define f_initialize_configuration     FORTRAN_NAME_(i_c,I_C)
#define f_initialize_geometry          FORTRAN_NAME_(i_g,I_G)
extern "C" {
   void f_initialize_configuration (const DOUBLE*,const INTEGER&);
   void f_initialize_geometry (const INTEGER&,const DOUBLE*);
}

#ifndef MIN
#  define MIN(a,b)                     (((a)<(b))?(a):(b))
#endif /* MIN */
#ifndef MAX
#  define MAX(a,b)                     (((a)>(b))?(a):(b))
#endif /* MAX */

class FlaggingSpecific :
public AMRFlagging<VectorType,FixupType,FlagType,DIM> {
   typedef AMRFlagging<VectorType,FixupType,FlagType,DIM> base;

   public:
      FlaggingSpecific(base::solver_type& solver) : base(solver) {
         base::AddCriterion(new ByValue<VectorType,FlagType,DIM>());
         base::AddCriterion(new ScaledGradient<VectorType,FlagType,DIM>());
         base::AddCriterion(new LimiterType<VectorType,FlagType,DIM>());
         base::AddCriterion(
            new AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver));
         base::AddCriterion(
            new RelativeError<VectorType,FixupType,FlagType,DIM>(solver));
#  ifdef f_flgout
         base::AddCriterion(new F77ByValue<VectorType,FlagType,DIM>(f_flgout));
         base::AddCriterion(
            new F77ScaledGradient<VectorType,FlagType,DIM>(f_flgout));
         base::AddCriterion(
            new F77LimiterType<VectorType,FlagType,DIM>(f_flgout));
         base::AddCriterion(
            new F77AbsoluteError<VectorType,FixupType,FlagType,DIM>(solver,
                                                                    f_flgout));
         base::AddCriterion(
            new F77RelativeError<VectorType,FixupType,FlagType,DIM>(solver,
                                                                    f_flgout));
#  endif /* f_flgout */
      }

      virtual void register_at(ControlDevice& Ctrl,const std::string& prefix) {
         std::string variable;
         int d,p;
         
         // Call the parent's register_at(ControlDevice&,const std::string&)
         // method
         base::register_at(Ctrl,prefix);

         for (p=0; p<NPATCHES; p++) {
            std::ostringstream patch;
            patch << p+1;
            for (d=0; d<DIM; d++) {
               std::ostringstream dimension;
               dimension << d+1;
               variable="GCellsMin" + patch.str() + "(" + dimension.str() + ")";
               RegisterAt(base::LocCtrl,variable.c_str(),_GCellsMin[p*DIM+d]);
               variable="GCellsMax" + patch.str() + "(" + dimension.str() + ")";
               RegisterAt(base::LocCtrl,variable.c_str(),_GCellsMax[p*DIM+d]);
            }
         }
      }

      virtual void register_at(ControlDevice& Ctrl) {
         register_at(Ctrl,"");
      }

      virtual void SetFlags(const int Time,const int Level,double t,double dt) {
         int p;

         // Call the parent's SetFlags(const int,const int,double,double) method
         base::SetFlags(Time,Level,t,dt);

         for (p=0; p<NPATCHES; p++) {
            Coords lbound(DIM,&(_GCellsMin[p*DIM])),step(DIM,1),
                   ubound(DIM,&(_GCellsMax[p*DIM]));

            BBox boundingBox(lbound*StepSize(base::GH(),0),
                             ubound*StepSize(base::GH(),0),
                             step*StepSize(base::GH(),0));
            boundingBox.refine(
               StepSize(base::GH(),0)/StepSize(base::GH(),Level));

            // FIXME: What does c stand for?
            forall (base::Flags(),Time,Level,c)
               base::Flags()(Time,Level,c).equals(BadPoint,boundingBox);
            end_forall
         }
      }

      ~FlaggingSpecific() {
         DeleteAllCriterions();
      }

   protected:
      int _GCellsMin[NPATCHES*DIM],_GCellsMax[NPATCHES*DIM];
};

class SolverSpecific :
public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
   typedef VectorType::InternalDataType DataType;
   typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
   typedef F77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;

   public:
      SolverSpecific(IntegratorSpecific& integ,
         base::initial_condition_type& init,
         base::boundary_conditions_type& bc) : base(integ, init, bc) {
         SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong,
                                                               f_restrict));
#  ifdef f_flgout
         SetFileOutput(
            new F77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,
                                                                    f_flgout));
#  else
         SetFileOutput(
            new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this));
#  endif /* f_flgout */
         SetFixup(new FixupSpecific());
         SetFlagging(new FlaggingSpecific(*this));
         AddGFM(new GhostFluidMethod<VectorType,DIM>(
            new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
            new F77GFMLevelSet<DataType,DIM>(f_lset)));
      }

      virtual void register_at(ControlDevice& Ctrl,const std::string& prefix) {
         ControlDevice device;

         // Call the parent's register_at(ControlDevice&,const std::string&)
         // method
         base::register_at(Ctrl,prefix);

         // Create a new ControlDevice for configuration options that are
         // specific to the geometry
         device=Ctrl.getSubDevice(prefix+"Geometry");
         // Register a new variable with the newly created ControlDevice
         RegisterAt(device,"TipRadius",_tipRadius);
      }

      virtual void register_at(ControlDevice& Ctrl) {
         // Call the overloaded register_at(ControlDevice&,const std::string&)
         // method
         register_at(Ctrl,"");
      }

      virtual void SetupData() {
         DOUBLE h=FLT_MAX;
         INTEGER refinement=1;
         int i;

         for (i=0; i<MaxLev; i++)
            refinement *= _RefineFactor[i];

         for (i=0; i<DIM; i++)
            h=MIN(h,(geom[2*i+1]-geom[2*i])/
                                      static_cast<DOUBLE>(refinement*shape[i]));
#  define NOPTIONS              2
         INTEGER n=NOPTIONS;
         DOUBLE options[NOPTIONS]={
                                    static_cast<DOUBLE>(_boundaryLayer),
                                    _tipRadius
                                    // _tipRadius*h
                                  };
#  undef NOPTIONS

         // Call the parent's SetupData() method
         base::SetupData();
         f_initialize_configuration (options,n);
         f_initialize_geometry (DIM,geom);
      }

      ~SolverSpecific() {
         DeleteGFM(_GFM[0]);
         delete _Flagging;
         delete _Fixup;
         delete _FileOutput;
      }
   protected:
      DOUBLE _tipRadius;
      INTEGER _boundaryLayer;