00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 2
00007
00008 #include "LBMProblem.h"
00009 #include "LBMD2Q9.h"
00010
00011 typedef double DataType;
00012 typedef LBMD2Q9<DataType> LBMType;
00013
00014 #define OWN_AMRSOLVER
00015 #define OWN_BOUNDARYCONDITION
00016 #include "LBMStdProblem.h"
00017 #include "LBMD2Q9SupplementalBC.h"
00018 #include "AMRGFMSolver.h"
00019 #include "LBMGFMBoundary.h"
00020 #include "AMRGFMInterpolation.h"
00021 #include "Interfaces/SchemeGFMFileOutput.h"
00022 #include "Interfaces/SchemeBrepOutput.h"
00023 #include "MultiStatCPTLevelSet.h"
00024
00025
00026 class BoundaryConditionsSpecific : public BoundaryConditionsSupplemental {
00027 public:
00028 BoundaryConditionsSpecific(LBMType &lbm) : BoundaryConditionsSupplemental(lbm) {}
00029 };
00030
00031
00032 class SolverSpecific :
00033 public AMRGFMSolver<MicroType,FixupType,FlagType,DIM> {
00034 typedef AMRGFMSolver<MicroType,FixupType,FlagType,DIM> base;
00035 typedef AMRGFMInterpolation<MicroType,FixupType,FlagType,DIM> interpolation_type;
00036 typedef MicroType::InternalDataType DataType;
00037 typedef SchemeBrepOutput<LBMType,DIM> brep_output_type;
00038 typedef SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM> output_type;
00039 typedef brep_output_type::geom_point_type geom_point_type;
00040 typedef MultiStatCPTLevelSet<DataType,DIM> cpt_type;
00041 typedef interpolation_type::point_type point_type;
00042 public:
00043 SolverSpecific(IntegratorSpecific& integ,
00044 base::initial_condition_type& init,
00045 base::boundary_conditions_type& bc) : base(integ, init, bc) {
00046 SetLevelTransfer(new LBMF77LevelTransfer<LBMType,DIM>(integ.Scheme(), f_prolong, f_restrict));
00047 SetFileOutput(new SchemeGFMFileOutput<LBMType,FixupType,FlagType,DIM>(*this,integ.Scheme()));
00048 SetFixup(new FixupSpecific(integ.Scheme()));
00049 SetFlagging(new FlaggingSpecific(*this,integ.Scheme()));
00050 AddGFM(new GhostFluidMethod<MicroType,DIM>(new LBMGFMBoundary<LBMType,DIM>(integ.LBM()),
00051 new cpt_type()));
00052 _BrepOutput = new brep_output_type(new interpolation_type(*this),integ.Scheme());
00053 IntegrateEvery = 1;
00054 moment_origin = 0.;
00055 }
00056
00057 ~SolverSpecific() {
00058 DeleteGFM(_GFM[0]);
00059 delete _LevelTransfer;
00060 delete _Flagging;
00061 delete _Fixup;
00062 delete _FileOutput;
00063 delete _BrepOutput;
00064 }
00065
00066 virtual void SetupData() {
00067 base::SetupData();
00068 base::NAMRTimeSteps = 1;
00069 base::AdaptBndTimeInterpolate = 0;
00070 base::Step[0].LastTime = 1.e37;
00071 base::Step[0].VariableTimeStepping = -1;
00072 base::Step[0].dtv[0] = ((LBMIntegrator<LBMType,DIM> &)Integrator_()).LBM().TimeScale();
00073 if (_BrepOutput) _BrepOutput->SetupData(base::PGH(), base::NGhosts());
00074 }
00075
00076 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00077 base::register_at(Ctrl,prefix);
00078 if (_BrepOutput) _BrepOutput->register_at(base::LocCtrl,prefix);
00079 RegisterAt(base::LocCtrl,"OutputDragEvery",IntegrateEvery);
00080 char VariableName[32];
00081 for (int d=0; d<base::Dim(); d++) {
00082 std::sprintf(VariableName,"MOrigin(%d)",d+1);
00083 RegisterAt(base::LocCtrl,VariableName,moment_origin[d]);
00084 }
00085 }
00086 virtual void register_at(ControlDevice& Ctrl) {
00087 base::register_at(Ctrl);
00088 }
00089
00090 virtual void update() {
00091 base::update();
00092 if (_BrepOutput) _BrepOutput->update();
00093 }
00094
00095 virtual void Advance(double& t, double& dt) {
00096 base::Advance(t,dt);
00097 if (StepsTaken(base::GH(),0)%IntegrateEvery == 0) {
00098 geom_point_type viscous_force(0.), viscous_moment(0.), total_force(0.), total_moment(0.);
00099
00100 int num_pos=0, num_elem=0;
00101 const DataType* vertices=0;
00102 const int* connections=0;
00103 if (!((cpt_type *)base::GFM(0).GetLevelSetP())->GetBrep(num_pos,vertices,num_elem,connections))
00104 return;
00105
00106 geom_point_type* centroids = new geom_point_type[num_elem];
00107 _BrepOutput->ElementCentroids(num_pos,vertices,num_elem,connections,centroids);
00108 point_type* xc = reinterpret_cast<point_type*>(centroids);
00109
00110
00111 DataType* dev = new DataType[3*num_elem];
00112 for (register int cnt=16; cnt<=18; cnt++) {
00113 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00114 int Time = CurrentTime(base::GH(),l);
00115 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00116 cnt, base::t[l]);
00117 }
00118
00119 _BrepOutput->Interpolation().PointsValues(base::Work(),base::t[0],num_elem,xc,dev+(cnt-16)*num_elem);
00120 }
00121 _BrepOutput->Interpolation().ArrayCombine(VizServer,3*num_elem,dev);
00122
00123 if (MY_PROC==VizServer) {
00124 geom_point_type* dev_tensor = new geom_point_type[2*num_elem];
00125 for (register int i=0; i<num_elem; i++) {
00126 dev_tensor[2*i ][0]=dev[ i]; dev_tensor[2*i ][1]=dev[ num_elem+i];
00127 dev_tensor[2*i+1][0]=dev[num_elem+i]; dev_tensor[2*i+1][1]=dev[2*num_elem+i];
00128 }
00129
00130 viscous_force = -_BrepOutput->IntegrateTensorNormal(dev_tensor,num_pos,vertices,num_elem,connections);
00131 viscous_moment = -_BrepOutput->IntegrateTensorMoment(dev_tensor,num_pos,vertices,num_elem,connections,moment_origin);
00132
00133 delete [] dev_tensor;
00134 }
00135 delete [] dev;
00136
00137
00138 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00139 int Time = CurrentTime(base::GH(),l);
00140 int press_cnt = base::Dim()+4;
00141 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00142 press_cnt, base::t[l]);
00143 }
00144
00145 geom_point_type pressure_force = _BrepOutput->IntegrateNormal(base::Work(),num_pos,vertices,num_elem,connections);
00146 geom_point_type pressure_moment = _BrepOutput->IntegrateMoment(base::Work(),num_pos,vertices,num_elem,connections,moment_origin);
00147
00148
00149 DataType* sig = new DataType[3*num_elem];
00150 for (register int cnt=23; cnt<=25; cnt++) {
00151 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00152 int Time = CurrentTime(base::GH(),l);
00153 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00154 cnt, base::t[l]);
00155 }
00156
00157 _BrepOutput->Interpolation().PointsValues(base::Work(),base::t[0],num_elem,xc,sig+(cnt-23)*num_elem);
00158 }
00159 _BrepOutput->Interpolation().ArrayCombine(VizServer,3*num_elem,sig);
00160
00161 if (MY_PROC==VizServer) {
00162 geom_point_type* sig_tensor = new geom_point_type[2*num_elem];
00163 for (register int i=0; i<num_elem; i++) {
00164 sig_tensor[2*i ][0]=sig[ i]; sig_tensor[2*i ][1]=sig[ num_elem+i];
00165 sig_tensor[2*i+1][0]=sig[num_elem+i]; sig_tensor[2*i+1][1]=sig[2*num_elem+i];
00166 }
00167
00168 total_force = -_BrepOutput->IntegrateTensorNormal(sig_tensor,num_pos,vertices,num_elem,connections);
00169 total_moment = -_BrepOutput->IntegrateTensorMoment(sig_tensor,num_pos,vertices,num_elem,connections,moment_origin);
00170
00171 delete [] sig_tensor;
00172 }
00173 delete [] sig;
00174 delete [] centroids;
00175
00176 if (MY_PROC==VizServer) {
00177 std::ofstream outfile;
00178 std::ostream* out;
00179 std::string fname = "drag.txt";
00180 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00181 out = new std::ostream(outfile.rdbuf());
00182 *out << base::t[0] << " " << total_force << " " << pressure_force << " " << viscous_force
00183 << " " << pressure_force+viscous_force << " " << total_moment << " " << pressure_moment << " " << viscous_moment
00184 << " " << pressure_moment+viscous_moment << std::endl;
00185 outfile.close();
00186 delete out;
00187 }
00188 }
00189 }
00190
00191 virtual void Output() {
00192 if (CurrentTime(base::GH(),0) == base::LastOutputTime()) return;
00193
00194 int num_vertices=0, num_edges=0;
00195 const DataType* vertices=0;
00196 const int* edges=0;
00197
00198 if (!((cpt_type *)base::GFM(0).GetLevelSetP())->GetBrep(num_vertices,vertices,num_edges,edges))
00199 return;
00200
00201 _BrepOutput->WriteOut(base::U(),num_vertices,vertices,num_edges,edges);
00202
00203 base::Output();
00204 }
00205
00206 protected:
00207 brep_output_type* _BrepOutput;
00208 int IntegrateEvery;
00209 geom_point_type moment_origin;
00210 };
00211
00212 #endif