00001
00002
00003 #ifndef AMROC_PROBLEM_H
00004 #define AMROC_PROBLEM_H
00005
00006 #define DIM 3
00007
00008 #include "LBMProblem.h"
00009 #include "LBMD3Q19.h"
00010
00011 typedef double DataType;
00012 typedef LBMD3Q19<DataType> LBMType;
00013
00014 #define OWN_AMRSOLVER
00015 #define OWN_BOUNDARYCONDITION
00016 #include "LBMStdProblem.h"
00017 #include "LBMD3Q19SupplementalBC.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[6*num_elem];
00112 for (register int cnt=19; cnt<=24; 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-19)*num_elem);
00120 }
00121 _BrepOutput->Interpolation().ArrayCombine(VizServer,6*num_elem,dev);
00122
00123 if (MY_PROC==VizServer) {
00124 geom_point_type* dev_tensor = new geom_point_type[3*num_elem];
00125 for (register int i=0; i<num_elem; i++) {
00126 dev_tensor[3*i ][0]=dev[ i]; dev_tensor[3*i ][1]=dev[ num_elem+i]; dev_tensor[3*i ][2]=dev[3*num_elem+i];
00127 dev_tensor[3*i+1][0]=dev[ num_elem+i]; dev_tensor[3*i+1][1]=dev[2*num_elem+i]; dev_tensor[3*i+1][2]=dev[4*num_elem+i];
00128 dev_tensor[3*i+2][0]=dev[3*num_elem+i]; dev_tensor[3*i+2][1]=dev[4*num_elem+i]; dev_tensor[3*i+2][2]=dev[5*num_elem+i];
00129 }
00130
00131 viscous_force = -_BrepOutput->IntegrateTensorNormal(dev_tensor,num_pos,vertices,num_elem,connections);
00132 viscous_moment = -_BrepOutput->IntegrateTensorMoment(dev_tensor,num_pos,vertices,num_elem,connections,moment_origin);
00133
00134 delete [] dev_tensor;
00135 }
00136 delete [] dev;
00137
00138
00139 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00140 int Time = CurrentTime(base::GH(),l);
00141 int press_cnt = base::Dim()+4;
00142 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00143 press_cnt, base::t[l]);
00144 }
00145
00146 geom_point_type pressure_force = _BrepOutput->IntegrateNormal(base::Work(),num_pos,vertices,num_elem,connections);
00147 geom_point_type pressure_moment = _BrepOutput->IntegrateMoment(base::Work(),num_pos,vertices,num_elem,connections,moment_origin);
00148
00149
00150 DataType* sig = new DataType[6*num_elem];
00151 for (register int cnt=28; cnt<=33; cnt++) {
00152 for (register int l=0; l<=FineLevel(base::GH()); l++) {
00153 int Time = CurrentTime(base::GH(),l);
00154 ((output_type*) _FileOutput)->Transform(base::U(), base::Work(), Time, l,
00155 cnt, base::t[l]);
00156 }
00157
00158 _BrepOutput->Interpolation().PointsValues(base::Work(),base::t[0],num_elem,xc,sig+(cnt-28)*num_elem);
00159 }
00160 _BrepOutput->Interpolation().ArrayCombine(VizServer,6*num_elem,sig);
00161
00162 if (MY_PROC==VizServer) {
00163 geom_point_type* sig_tensor = new geom_point_type[3*num_elem];
00164 for (register int i=0; i<num_elem; i++) {
00165 sig_tensor[3*i ][0]=sig[ i]; sig_tensor[3*i ][1]=sig[ num_elem+i]; sig_tensor[3*i ][2]=sig[3*num_elem+i];
00166 sig_tensor[3*i+1][0]=sig[ num_elem+i]; sig_tensor[3*i+1][1]=sig[2*num_elem+i]; sig_tensor[3*i+1][2]=sig[4*num_elem+i];
00167 sig_tensor[3*i+2][0]=sig[3*num_elem+i]; sig_tensor[3*i+2][1]=sig[4*num_elem+i]; sig_tensor[3*i+2][2]=sig[5*num_elem+i];
00168 }
00169
00170 total_force = -_BrepOutput->IntegrateTensorNormal(sig_tensor,num_pos,vertices,num_elem,connections);
00171 total_moment = -_BrepOutput->IntegrateTensorMoment(sig_tensor,num_pos,vertices,num_elem,connections,moment_origin);
00172
00173 delete [] sig_tensor;
00174 }
00175 delete [] sig;
00176 delete [] centroids;
00177
00178 if (MY_PROC==VizServer) {
00179 std::ofstream outfile;
00180 std::ostream* out;
00181 std::string fname = "drag.txt";
00182 outfile.open(fname.c_str(), std::ios::out | std::ios::app);
00183 out = new std::ostream(outfile.rdbuf());
00184 *out << base::t[0] << " " << total_force << " " << pressure_force << " " << viscous_force
00185 << " " << pressure_force+viscous_force << " " << total_moment << " " << pressure_moment << " " << viscous_moment
00186 << " " << pressure_moment+viscous_moment << std::endl;
00187 outfile.close();
00188 delete out;
00189 }
00190 }
00191 }
00192
00193 virtual void Output() {
00194 if (CurrentTime(base::GH(),0) == base::LastOutputTime()) return;
00195
00196 int num_vertices=0, num_edges=0;
00197 const DataType* vertices=0;
00198 const int* edges=0;
00199
00200 if (!((cpt_type *)base::GFM(0).GetLevelSetP())->GetBrep(num_vertices,vertices,num_edges,edges))
00201 return;
00202
00203 _BrepOutput->WriteOut(base::U(),num_vertices,vertices,num_edges,edges);
00204
00205 base::Output();
00206 }
00207
00208 protected:
00209 brep_output_type* _BrepOutput;
00210 int IntegrateEvery;
00211 geom_point_type moment_origin;
00212 };
00213
00214 #endif