diff --git a/demo/src/navierStokes.h b/demo/src/navierStokes.h
index 771f050117fdf6441c3da99b6cfe770659d5d7ec..ed482349b19a9c55a8f29e694b45e0507863cb08 100644
--- a/demo/src/navierStokes.h
+++ b/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;
diff --git a/demo/src/navierStokes_diffuseDomain.cc b/demo/src/navierStokes_diffuseDomain.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e80d05bf5c116313e86790f0133963be472c2ef2
--- /dev/null
+++ b/demo/src/navierStokes_diffuseDomain.cc
@@ -0,0 +1,196 @@
+#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;
+};
diff --git a/demo/src/pfc.cc b/demo/src/pfc.cc
index 387de288792577d9dc65dd0e1946eda33ce50be1..74cfeaeea21bcc10b71f866a154d07041bfb7aeb 100644
--- a/demo/src/pfc.cc
+++ b/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);
+// }
+// }
}
};