navierStokes dome

demo/src/navierStokes.h

@@ -39,7 +39,39 @@ public:
     for (size_t i = 0; i < vertices.size()-1; i++)
       result = std::min(result, meshconv::distance_point_line_2d(x.begin(), vertices[i].begin(), vertices[i+1].begin()));
     return result * (meshconv::point_in_polygon(x.begin(), vertices) ? -1.0 : 1.0);
-  };
+  }
+  
+  void refine(int np)
+  {
+    std::vector<WorldVector<double> > newVertices;
+    
+    for (size_t i = 0; i < vertices.size()-1; i++) {
+      for (size_t j = 0; j < np-1; j++) {
+	double lambda = static_cast<double>(j)/static_cast<double>(np - 1.0);
+	WorldVector<double> p = lambda*vertices[i+1] + (1.0-lambda)*vertices[i];
+	newVertices.push_back(p);
+      }
+    }
+    swap(vertices, newVertices);
+  }
+  
+  void move(const DOFVector<WorldVector<double> >* velocity)
+  {
+    for (size_t i = 0; i < vertices.size()-1; i++) {
+      WorldVector<double> shift = evalAtPoint(*velocity, vertices[i]);
+      vertices[i] += shift;
+    }
+    vertices[vertices.size()-1] = vertices[0];
+  }
+  
+  void move(AbstractFunction<WorldVector<double>, WorldVector<double> >* velocity)
+  {
+    for (size_t i = 0; i < vertices.size()-1; i++) {
+      WorldVector<double> shift = (*velocity)(vertices[i]);
+      vertices[i] += shift;
+    }
+    vertices[vertices.size()-1] = vertices[0];
+  }
 
 private:
   std::vector<WorldVector<double> > vertices;

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 "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_, double tau_) : 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));
+  }
+  
+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, double* time, double tau, WorldVector<double> fixedPoint_)
+  : 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;
+  }
+  
+private:
+  AbstractFunction<double, double>* fct;
+  WorldVector<double> fixedPoint;
+};
+
+class NS_Channel : public NavierStokesPhase_TaylorHood
+{
+public:
+  typedef NavierStokesPhase_TaylorHood super;
+  
+public:
+  NS_Channel(std::string name_) : super(name_) {}
+  
+  ~NS_Channel()
+  {
+    delete phaseField;
+  }
+  
+  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));
+    
+    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, adaptInfo->getTimePtr(), adaptInfo->getTimestep(), fixedPoint);
+    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);
+
+    WorldVector<double> nullVec; nullVec.set(0.0);
+    AbstractFunction<WorldVector<double>, WorldVector<double> > *zeroVec = new AMDiS::Const<WorldVector<double>, WorldVector<double> >(nullVec);
+    super::setBcFct(zeroVec);
+    
+    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;
+  AbstractFunction<WorldVector<double>, WorldVector<double> >* 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;
+};

demo/src/pfc.cc

@@ -44,17 +44,22 @@ struct PFC_Demo : public PhaseFieldCrystal
 
     // random data
     prob->getSolution()->getDOFVector(0)->interpol(new RandomValues(density, amplitude));
+    
+    
+    std::cout << "Number of degrees of freedom: "
+	      << prob->getSolution()->getDOFVector(0)->getUsedSize()
+	      << std::endl;
   }
 
   // add periodic boundary condition on all boundaries
   void fillBoundaryConditions()
   {
-    for (size_t i = 0; i < prob->getNumComponents(); i++) {
-    for (size_t j = 0; j < prob->getNumComponents(); j++) {
-      prob->addPeriodicBC(-1, i, j);
-      prob->addPeriodicBC(-2, i, j);
-    }
-    }
+//     for (size_t i = 0; i < prob->getNumComponents(); i++) {
+//     for (size_t j = 0; j < prob->getNumComponents(); j++) {
+//       prob->addPeriodicBC(-1, i, j);
+//       prob->addPeriodicBC(-2, i, j);
+//     }
+//     }
   }
 };