2 navierstokes demos added

extensions/demo/src/navierStokes_diffuseDomain.cc

+#include "AMDiS.h"
+#include "NavierStokesPhase_TaylorHood.h"
+#include "navierStokes.h"
+// #include "time/ExtendedRosenbrockAdaptInstationary.h"
+#include "Refinement.h"
+#include "MeshFunction_Level.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+class NS_Channel : public NavierStokesPhase_TaylorHood
+{
+public:
+  typedef NavierStokesPhase_TaylorHood super;
+  
+public:
+  NS_Channel(std::string name_) : super(name_), phaseField(NULL), inflowBC(NULL) {}
+  
+  ~NS_Channel()
+  {
+    if (phaseField != NULL) {
+      delete phaseField;
+      phaseField = NULL;
+    }
+    if (inflowBC != NULL) {
+      delete inflowBC;
+      inflowBC = NULL;
+    }
+  }
+  
+  void initData()
+  { FUNCNAME("NS_Channel::initData()");
+  
+    super::initData();
+  
+    phaseField = new DOFVector<double>(getFeSpace(0), "phaseField"); 
+    phaseField->set(1.0);
+    super::setPhase(phaseField);   
+    
+    double H = 4.1;
+    double Um = 1.5;
+    Parameters::get("mesh->H",H);
+    Parameters::get("ns->Um",Um);
+    inflowBC = new InflowBC(H,Um);
+  }
+
+  void solveInitialProblem(AdaptInfo *adaptInfo)
+  { FUNCNAME("NS_Channel::solveInitialProblem()");
+  
+    super::solveInitialProblem(adaptInfo);
+
+    double radius = 0.1;
+    WorldVector<double> center;
+    Parameters::get("obstacle->center",center);
+    Parameters::get("obstacle->radius",radius);
+    obstacle = new Circle(radius, center);
+
+    refFunction = new SignedDistRefinement(getMesh());
+    refinement = new RefinementLevelCoords2(getFeSpace(), refFunction, obstacle);
+
+    // initial refinement
+    refinement->refine(10);    
+    phaseField->interpol(new SignedDistFctToPhaseField(getEpsilon(), obstacle, -3.0));
+    
+    if (inflowBC)
+      inflowBC->setTimePtr(adaptInfo->getTimePtr());
+  }
+  
+protected:
+  
+  void fillBoundaryConditions()
+  { FUNCNAME("NS_Channel::fillBoundaryConditions()");
+    
+    WorldVector<double> zeroVec; zeroVec.set(0.0);
+    AbstractFunction<WorldVector<double>, WorldVector<double> > *zeroVecFct = new AMDiS::Const<WorldVector<double>, WorldVector<double> >(zeroVec);
+
+    super::setBcFct(zeroVecFct);    
+    super::fillBoundaryConditions();
+    
+    
+    // +------ 5 ------+
+    // |               |
+    // 2   # <--1      3
+    // |               |
+    // +------ 4 ------+
+    
+    AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
+    
+    /// at rigid wall: no-slip boundary condition
+    for (size_t i = 0; i < dow; i++) {
+      getProblem(0)->addDirichletBC(4, i, i, zero);
+      getProblem(0)->addDirichletBC(5, i, i, zero);
+    }
+    /// dirichlet bc for pressure at one DOF
+//     getProblem(0)->addSingularDirichletBC(0, dow, dow, *zero);
+
+
+    /// at left wall: prescribed velocity    
+    if (inflowBC) {
+      getProblem(0)->addDirichletBC(2, 0, 0, inflowBC);
+      getProblem(0)->addDirichletBC(2, 1, 1, zero);
+    }
+  }
+  
+private:
+  
+  DOFVector<double>* phaseField;
+  AbstractFunction<double, WorldVector<double> >* obstacle;
+  
+  SignedDistRefinement* refFunction;
+  RefinementLevelCoords2* refinement;
+  
+  InflowBC* inflowBC;
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+  AMDiS::init(argc, argv);
+
+  NS_Channel nsProb("ns");
+  nsProb.initialize(INIT_ALL | INIT_EXACT_SOLUTION);
+
+  // Adapt-Infos
+  AdaptInfo adaptInfo("adapt", nsProb.getNumComponents());
+
+  // adaption loop - solve ch-prob and ns-prob
+  AdaptInstationary adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
+//   ExtendedRosenbrockAdaptInstationary<NS_Channel> adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
+  
+  // scale Mesh
+  WorldVector<double> scale; scale[0] = 25.0; scale[1] = 4.1;
+  Helpers::scaleMesh(nsProb.getMesh(), scale);
+
+  ptime start_time = microsec_clock::local_time();
+  nsProb.initTimeInterface();
+  int error_code = adaptInstat.adapt(); 
+  time_duration td = microsec_clock::local_time()-start_time;
+
+  MSG("elapsed time= %d sec\n", td.total_seconds());
+
+  AMDiS::finalize();
+
+  return error_code;
+};

extensions/demo/src/navierStokes_diffuseDomain2.cc

+#include "AMDiS.h"
+#include "NavierStokesPhase_TaylorHood.h"
+#include "navierStokes.h"
+// #include "time/ExtendedRosenbrockAdaptInstationary.h"
+#include "Refinement.h"
+#include "MeshFunction_Level.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+
+#include "boost/date_time/posix_time/posix_time.hpp"
+
+using namespace AMDiS;
+using namespace boost::posix_time;
+
+struct BeamDisplacement1d : AbstractFunction<double, double>
+{
+  BeamDisplacement1d(double P_, double L_, double EI_, double* time_ = NULL, double* tau_ = NULL) : P(P_), L(L_), EI(EI_), time(time_), tau(tau_) {}
+
+  double operator()(const double& x) const
+  {
+    return (*time < DBL_TOL ? 0.0 : (sin(*time)-sin(*time - *tau))*P*sqr(x)*(3.0*L-x)/(6.0*EI));
+  }
+
+  void setTimePtr(double* time_)
+  {
+    time = time_;
+  }
+
+  void setTauPtr(double* tau_)
+  {
+    tau = tau_;
+  }
+  
+private:
+  double P;
+  double L;
+  double EI;
+  double* time;
+  double* tau;
+};
+
+struct BeamDisplacement : AbstractFunction<WorldVector<double>, WorldVector<double> >
+{
+  BeamDisplacement(double P, double L, double EI, WorldVector<double> fixedPoint_, double* time = NULL, double* tau = NULL)
+  : fct(new BeamDisplacement1d(P, L, EI, time, tau)), fixedPoint(fixedPoint_) {}
+  
+  WorldVector<double> operator()(const WorldVector<double>& x) const
+  {
+    WorldVector<double> displacement; displacement[0] = 0.0;
+    displacement[1] = (x[0] >= fixedPoint[0] ? (*fct)(x[0] - fixedPoint[0]) : 0.0);
+    return displacement;
+  }
+
+  void setTimePtr(double* time_)
+  {
+    fct->setTimePtr(time_);
+  }
+
+  void setTauPtr(double* tau_)
+  {
+    fct->setTauPtr(tau_);
+  }
+
+private:
+  BeamDisplacement1d* fct;
+  WorldVector<double> fixedPoint;
+};
+
+class NS_Channel : public NavierStokesPhase_TaylorHood
+{
+public:
+  typedef NavierStokesPhase_TaylorHood super;
+  
+public:
+  NS_Channel(std::string name_) : super(name_), phaseField(NULL) {}
+  
+  ~NS_Channel()
+  {
+    if (phaseField != NULL) {
+      delete phaseField;
+      phaseField = NULL;
+    }
+  }
+  
+  void initData()
+  { FUNCNAME("NS_Channel::initData()");
+  
+    super::initData();
+  
+    phaseField = new DOFVector<double>(getFeSpace(0), "phaseField"); 
+    phaseField->set(1.0);
+    super::setPhase(phaseField);   
+  }
+
+  void solveInitialProblem(AdaptInfo *adaptInfo)
+  { FUNCNAME("NS_Channel::solveInitialProblem()");
+  
+    super::solveInitialProblem(adaptInfo);
+
+    size_t nVertices = 0;
+    WorldVector<double> x;
+    Parameters::get("obstacle->num vertices",nVertices);
+    std::vector<WorldVector<double> > v(nVertices,x);
+    for (size_t i = 0; i < nVertices; i++)
+      Parameters::get("obstacle->vertex["+boost::lexical_cast<std::string>(i)+"]",v[i]);
+    v.push_back(v[0]);
+    
+    obstacle = new Polygon(v);
+
+    refFunction = new SignedDistRefinement(getMesh());
+    refinement = new RefinementLevelCoords2(
+      getFeSpace(),
+      refFunction,
+      obstacle);
+
+    // initial refinement
+    refinement->refine(10);
+    
+    phaseField->interpol(new SignedDistFctToPhaseField(getEpsilon(), obstacle, -3.0));
+    beamDisplacement->setTimePtr(adaptInfo->getTimePtr());
+    beamDisplacement->setTauPtr(adaptInfo->getTimestepPtr());
+    obstacle->refine(10);
+  }
+  
+protected:
+  
+  void fillBoundaryConditions()
+  { FUNCNAME("NS_Channel::fillBoundaryConditions()");
+    
+    AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
+    size_t dow = Global::getGeo(WORLD);
+
+    double P = 0.1, L = 2.0, EI = 1.0;
+    WorldVector<double> fixedPoint; fixedPoint[0] = 2.5; fixedPoint[1] = 2.0;
+    Parameters::get("obstacle->P",P);
+    Parameters::get("obstacle->L",L);
+    Parameters::get("obstacle->EI",EI);
+    Parameters::get("obstacle->fixed point",fixedPoint);
+    beamDisplacement = new BeamDisplacement(P, L, EI, fixedPoint);
+    super::setBcFct(beamDisplacement);
+    
+    super::fillBoundaryConditions();
+    
+    
+    // +------ 5 ------+
+    // |               |
+    // 2   # <--1      3
+    // |               |
+    // +------ 4 ------+
+    
+    /// at rigid wall: no-slip boundary condition
+    for (size_t i = 0; i < dow; i++) {
+      getProblem(0)->addDirichletBC(4, i, i, zero);
+      getProblem(0)->addDirichletBC(5, i, i, zero);
+    }
+    /// dirichlet bc for pressure at one DOF
+//     getProblem(0)->addSingularDirichletBC(0, dow, dow, *zero);
+
+    double H = 4.1;
+    double Um = 1.5;
+    Parameters::get("mesh->H",H);
+    Parameters::get("ns->Um",Um);
+
+    /// at left wall: prescribed velocity
+//     getProblem(0)->addDirichletBC(2, 0, 0, new InflowBC(H,Um));
+//     getProblem(0)->addDirichletBC(2, 1, 1, zero);
+  }
+  
+  void initTimestep(AdaptInfo* adaptInfo)
+  {
+    super::initTimestep(adaptInfo);
+    
+    obstacle->move(beamDisplacement);
+    refinement->refine(5);
+    
+    phaseField->interpol(new SignedDistFctToPhaseField(getEpsilon(), obstacle, -3.0));
+  } 
+  
+private:
+  
+  DOFVector<double>* phaseField;
+  BeamDisplacement* beamDisplacement;
+  Polygon* obstacle;
+  
+  SignedDistRefinement* refFunction;
+  RefinementLevelCoords2* refinement;
+};
+
+
+int main(int argc, char** argv)
+{ FUNCNAME("main");
+
+  AMDiS::init(argc, argv);
+
+  NS_Channel nsProb("ns");
+  nsProb.initialize(INIT_ALL | INIT_EXACT_SOLUTION);
+
+  // Adapt-Infos
+  AdaptInfo adaptInfo("adapt", nsProb.getNumComponents());
+
+  // adaption loop - solve ch-prob and ns-prob
+  AdaptInstationary adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
+//   ExtendedRosenbrockAdaptInstationary<NS_Channel> adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
+  
+  // scale Mesh
+  WorldVector<double> scale; scale[0] = 25.0; scale[1] = 4.1;
+  Helpers::scaleMesh(nsProb.getMesh(), scale);
+
+  ptime start_time = microsec_clock::local_time();
+  nsProb.initTimeInterface();
+  int error_code = adaptInstat.adapt(); 
+  time_duration td = microsec_clock::local_time()-start_time;
+
+  MSG("elapsed time= %d sec\n", td.total_seconds());
+
+  AMDiS::finalize();
+
+  return error_code;
+};