00001
00002
00003
00004
00005
00006 #ifndef AMROC_PROBLEM_H
00007 #define AMROC_PROBLEM_H
00008
00009 #include <cmath>
00010 #define DIM 3
00011
00012 #include "WENOProblem.h"
00013
00014 #define f_combl_read FORTRAN_NAME(cblread, CBLREAD)
00015
00016 extern "C" {
00017 void f_combl_read(INTEGER& NSph, DOUBLE xm[], DOUBLE rd[], INTEGER Np[]);
00018 }
00019
00020 #define OWN_GFMAMRSOLVER
00021 #include "WENOStdGFMProblem.h"
00022 #include "AMRGFMInterpolation.h"
00023
00024 class SolverSpecific :
00025 public AMRGFMSolver<VectorType,FixupType,FlagType,DIM> {
00026 typedef VectorType::InternalDataType DataType;
00027 typedef AMRGFMSolver<VectorType,FixupType,FlagType,DIM> base;
00028 typedef AMRGFMInterpolation<VectorType,FixupType,FlagType,DIM> interpolation_type;
00029 typedef interpolation_type::point_type point_type;
00030 typedef WENOF77GFMFileOutput<VectorType,FixupType,FlagType,DIM> output_type;
00031 public:
00032 SolverSpecific(IntegratorSpecific& integ,
00033 base::initial_condition_type& init,
00034 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00035 SetLevelTransfer(new F77LevelTransfer<VectorType,DIM>(f_prolong, f_restrict));
00036 #ifdef f_flgout
00037 SetFileOutput(new WENOF77GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this,f_flgout,f_bounds));
00038 #else
00039 SetFileOutput(new GFMFileOutput<VectorType,FixupType,FlagType,DIM>(*this));
00040 #endif
00041 SetFixup(new FixupSpecific());
00042 SetFlagging(new FlaggingSpecific(*this));
00043 AddGFM(new GhostFluidMethod<VectorType,DIM>(
00044 new F77GFMBoundary<VectorType,DIM>(f_ibndrfl,f_itrans),
00045 new F77GFMLevelSet<DataType,DIM>(f_lset)));
00046 _Interpolation = new interpolation_type(*this);
00047 }
00048
00049 ~SolverSpecific() {
00050 DeleteGFM(_GFM[0]);
00051 delete _Flagging;
00052 delete _Fixup;
00053 delete _Interpolation;
00054 }
00055
00056 virtual void SetupData() {
00057 base::SetupData();
00058 _Interpolation->SetupData(base::PGH(), base::NGhosts());
00059 }
00060
00061
00062
00063
00064 virtual double Tick(int VariableTimeStepping, const double dtv[],
00065 const double cflv[], int& Rejections) {
00066
00067 double cfl = base::Tick(VariableTimeStepping, dtv, cflv, Rejections);
00068
00069 int me = MY_PROC;
00070 double pi = 4.0*std::atan(1.0);
00071
00072
00073 int NSph, Np[5];
00074 double xm[15], rd[5];
00075 point_type sum[5];
00076 f_combl_read(NSph,xm,rd,Np);
00077
00078
00079 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00080 int Time = CurrentTime(base::GH(),l);
00081 int press_cnt = base::Dim()+4;
00082 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00083 press_cnt, base::t[l]);
00084 }
00085
00086 for (int ns=0; ns<NSph; ns++) {
00087 int Npoints = 8*Np[ns]*Np[ns];
00088 DataType* p = new DataType[Npoints];
00089 point_type* xc = new point_type[Npoints];
00090 point_type* normals = new point_type[Npoints];
00091 point_type* ds = new point_type[Npoints];
00092 DataType* dist = (DataType*) 0;
00093 double dalpha = 0.5*pi/Np[ns];
00094 double dbeta = 0.5*pi/Np[ns];
00095 register int n = 0;
00096
00097
00098 for (register int l=0; l<4*Np[ns]; l++) {
00099 double alpha = dalpha*(l+0.5);
00100 for (register int m=0; m<2*Np[ns]; m++) {
00101 double beta = dalpha*(m+0.5);
00102 xc[n](0) = xm[base::Dim()*ns ]+rd[ns]*std::cos(alpha)*std::sin(beta);
00103 xc[n](1) = xm[base::Dim()*ns+1]+rd[ns]*std::sin(alpha)*std::sin(beta);
00104 xc[n](2) = xm[base::Dim()*ns+2]+rd[ns]*std::cos(beta);
00105 ds[n] = pi*rd[ns]*rd[ns]*dbeta*std::sin(beta)/(2.*Np[ns]);
00106 n++;
00107 }
00108 }
00109
00110
00111 _Interpolation->PointsValues(base::Work(),base::t[0],Npoints,xc,p);
00112 _Interpolation->ArrayCombine(VizServer,Npoints,p);
00113
00114
00115 _Interpolation->PointsGeom(base::GFM(0).Phi(),base::t[0],Npoints,xc,normals,dist);
00116 _Interpolation->ArrayCombine(VizServer,base::Dim()*Npoints,normals[0].data());
00117
00118
00119 if (me == VizServer) {
00120 sum[ns] = 0.0;
00121 for (n=0; n<Npoints; n++)
00122 sum[ns] += p[n]*normals[n]*ds[n];
00123 }
00124
00125 delete [] p;
00126 delete [] xc;
00127 delete [] normals;
00128 delete [] ds;
00129 }
00130
00131
00132 if (me == VizServer) {
00133 std::ofstream outfile;
00134 std::ostream* out;
00135 std::string fname = "drag.dat";
00136 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00137 out = new std::ostream(outfile.rdbuf());
00138 *out << base::t[0];
00139 for (int ns=0; ns<NSph; ns++)
00140 for (int d=0; d<base::Dim(); d++)
00141 *out << " " << sum[ns](d);
00142 *out << std::endl;
00143 outfile.close();
00144 delete out;
00145 }
00146
00147 return cfl;
00148 }
00149
00150 protected:
00151 interpolation_type* _Interpolation;