diff --git a/extensions/ExtendedProblemStat.h b/extensions/ExtendedProblemStat.h
index 2b1c8e1231f6621f86bb72838f9bdeb1bfa6a0fc..06359ad6b709b27c2242ccb9dd8395016e3f3010 100644
--- a/extensions/ExtendedProblemStat.h
+++ b/extensions/ExtendedProblemStat.h
@@ -129,7 +129,7 @@ public:
if (asmMatrix && (singularDirichletBC.size() > 0 || manualPeriodicBC.size() > 0)) {
solverMatrix.setMatrix(*getSystemMatrix());
}
-
+// writeMatrix("matrix.mat");
#endif
}
diff --git a/extensions/POperators_ZOT.h b/extensions/POperators_ZOT.h
index 8a950b2af8d8eb898b2e1ba1c9000f4cc30e3029..5306f3635742bc6e957f797c8962bf78f38df502 100644
--- a/extensions/POperators_ZOT.h
+++ b/extensions/POperators_ZOT.h
@@ -22,9 +22,9 @@ void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad=NULL);
-inline void getC(const ElInfo *, int nPoints, ElementVector &C);
+void getC(const ElInfo *, int nPoints, ElementVector &C);
-inline void eval(int nPoints,
+void eval(int nPoints,
const mtl::dense_vector<double>&,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
@@ -76,16 +76,16 @@ void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad=NULL);
-inline void getC(const ElInfo *, int nPoints, ElementVector &C);
+void getC(const ElInfo *, int nPoints, ElementVector &C);
-inline void eval(int nPoints,
+void eval(int nPoints,
const mtl::dense_vector<double>&,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double opFactor);
-inline double f(const int iq) const;
+double f(const int iq) const;
void setFactor(const double fac_) { fac=fac_; }
protected:
@@ -102,26 +102,26 @@ WorldVector<double>* facVec;
class Pow3_ZOT : public ZeroOrderTerm
{
public:
- Pow3_ZOT(DOFVectorBase<double> *rhoDV_, double fac_=1.0);
- void initElement(const ElInfo* elInfo,
- SubAssembler* subAssembler,
- Quadrature *quad = NULL);
- void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
- SubAssembler*,
- Quadrature *quad = NULL);
-inline void getC(const ElInfo *, int nPoints, ElementVector &C);
+ Pow3_ZOT(DOFVectorBase<double> *rhoDV_, double fac_=1.0);
+ void initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+ void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+ SubAssembler*,
+ Quadrature *quad = NULL);
+ void getC(const ElInfo *, int nPoints, ElementVector &C);
-inline void eval(int nPoints,
- const mtl::dense_vector<double>&,
- const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
- const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
- mtl::dense_vector<double>& result,
- double opFactor);
- inline double f(const int iq) const;
+ void eval(int nPoints,
+ const mtl::dense_vector<double>&,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double opFactor);
+ double f(const int iq) const;
protected:
- DOFVectorBase<double> *rhoDV;
- mtl::dense_vector<double> rho;
- double fac;
+ DOFVectorBase<double> *rhoDV;
+ mtl::dense_vector<double> rho;
+ double fac;
};
/* -------------------------------------------------------------- */
@@ -137,7 +137,7 @@ public:
void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler*,
Quadrature *quad = NULL);
- inline double f(const int iq) const;
+ double f(const int iq) const;
protected:
DOFVectorBase<double> *rhoDV;
mtl::dense_vector<double> rho;
@@ -148,26 +148,26 @@ protected:
class Pow2_ZOT : public ZeroOrderTerm
{
public:
- Pow2_ZOT(DOFVectorBase<double> *rhoDV_, double fac_=1.0);
- void initElement(const ElInfo* elInfo,
- SubAssembler* subAssembler,
- Quadrature *quad = NULL);
- void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
- SubAssembler*,
- Quadrature *quad = NULL);
-inline void getC(const ElInfo *, int nPoints, ElementVector &C);
+ Pow2_ZOT(DOFVectorBase<double> *rhoDV_, double fac_=1.0);
+ void initElement(const ElInfo* elInfo,
+ SubAssembler* subAssembler,
+ Quadrature *quad = NULL);
+ void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
+ SubAssembler*,
+ Quadrature *quad = NULL);
+ void getC(const ElInfo *, int nPoints, ElementVector &C);
-inline void eval(int nPoints,
- const mtl::dense_vector<double>&,
- const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
- const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
- mtl::dense_vector<double>& result,
- double opFactor);
- inline double f(const int iq) const;
+ void eval(int nPoints,
+ const mtl::dense_vector<double>&,
+ const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
+ const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
+ mtl::dense_vector<double>& result,
+ double opFactor);
+ double f(const int iq) const;
protected:
- DOFVectorBase<double> *rhoDV;
- mtl::dense_vector<double> rho;
- double fac;
+ DOFVectorBase<double> *rhoDV;
+ mtl::dense_vector<double> rho;
+ double fac;
};
/* -------------------------------------------------------------- */
@@ -183,7 +183,7 @@ public:
void initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler*,
Quadrature *quad = NULL);
- inline double f(const int iq) const;
+ double f(const int iq) const;
protected:
DOFVectorBase<double> *rhoDV;
mtl::dense_vector<double> rho;
diff --git a/extensions/SignedDistFunctors.h b/extensions/SignedDistFunctors.h
index ae5e72d927f5805c84a110e6dfee41094456ba50..87bb1501b8414ec7256dc036ea9d83facc88f823 100644
--- a/extensions/SignedDistFunctors.h
+++ b/extensions/SignedDistFunctors.h
@@ -457,10 +457,15 @@ public:
double operator()(const WorldVector<double>& x) const
{
double result=0.0;
- PerturbedCircleRadius perturbedCircle(radius, 1.0, 3, rotation);
+ PerturbedCircleRadius perturbedCircle(radius, 0.7, 3, rotation);
result = signedDist2D(x, midPoint, &perturbedCircle, 1.e-6);
return -result;
- };
+ }
+
+ void setRotation(double rotation_)
+ {
+ rotation = rotation_;
+ }
private:
diff --git a/extensions/SingularDirichletBC.h b/extensions/SingularDirichletBC.h
index e23b2401f3755de5c3ebc4ebfbc1ef8db0b72534..cb9b28eabbe44c166f00c7ce15fb1fac8e15944e 100644
--- a/extensions/SingularDirichletBC.h
+++ b/extensions/SingularDirichletBC.h
@@ -185,88 +185,30 @@ namespace details {
* The 'nr' and 'meshIndicator' describe one boundary of the mesh and to each DegreeOfFreedom
* the 'periodicMap' associates a second DegreeOfFreedom.
**/
- inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
- BoundaryType nr,
- AbstractFunction<bool, WorldVector<double> >* meshIndicator,
- AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
- std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
- {
- std::set<DegreeOfFreedom> indices;
- details::getBoundaryIndices(feSpace, nr, meshIndicator, indices);
-
- association.clear();
-
- DOFVector<WorldVector<double> > coords(feSpace, "coords");
- feSpace->getMesh()->getDofIndexCoords(coords);
- std::set<DegreeOfFreedom>::iterator it;
- for (it = indices.begin(); it != indices.end(); it++)
- {
- DegreeOfFreedom idx2;
- WorldVector<double> p;
- p = (*periodicMap)(coords[*it]);
- TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
-
- association.push_back(std::make_pair(*it, idx2));
- }
- }
- inline bool getDofIdxAtPoint(const FiniteElemSpace* feSpace,
- WorldVector<double> &p,
- DegreeOfFreedom &idx,
- ElInfo *oldElInfo = NULL,
- bool useOldElInfo = false)
- {
- FUNCNAME("DOFVector::getDofIdxAtPoint()");
-
- Mesh *mesh = feSpace->getMesh();
- const BasisFunction *basFcts = feSpace->getBasisFcts();
-
- int dim = mesh->getDim();
- int numBasFcts = basFcts->getNumber();
-
- std::vector<DegreeOfFreedom> localIndices(numBasFcts);
- DimVec<double> lambda(dim, NO_INIT);
-
- ElInfo *elInfo = mesh->createNewElInfo();
- idx = 0;
- bool inside = false;
-
- if (oldElInfo && useOldElInfo && oldElInfo->getMacroElement()) {
- inside = mesh->findElInfoAtPoint(p, elInfo, lambda, oldElInfo->getMacroElement(), NULL, NULL);
- delete oldElInfo;
- } else {
- inside = mesh->findElInfoAtPoint(p, elInfo, lambda, NULL, NULL, NULL);
- }
-
- if (oldElInfo)
- oldElInfo = elInfo;
-
- if (!inside)
- return false;
-
-
- basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
-
-
- WorldVector<double> coord;
- int minIdx = -1;
- double minDist = 1.e15;
-
- for (int i = 0; i < numBasFcts; i++) {
- elInfo->coordToWorld(*(basFcts->getCoords(i)), coord);
- WorldVector<double> dist = coord - p;
- double newDist = norm(dist);
- if (newDist < minDist) {
- minIdx = i;
- minDist = newDist;
- }
- }
-
- if (minIdx >= 0)
- idx = localIndices[minIdx];
-
- if(!oldElInfo) delete elInfo;
- return inside;
- }
+// inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
+// BoundaryType nr,
+// AbstractFunction<bool, WorldVector<double> >* meshIndicator,
+// AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
+// std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
+// {
+// std::set<DegreeOfFreedom> indices;
+// details::getBoundaryIndices(feSpace, nr, meshIndicator, indices);
+//
+// association.clear();
+//
+// DOFVector<WorldVector<double> > coords(feSpace, "coords");
+// feSpace->getMesh()->getDofIndexCoords(coords);
+// std::set<DegreeOfFreedom>::iterator it;
+// for (it = indices.begin(); it != indices.end(); it++)
+// {
+// DegreeOfFreedom idx2;
+// WorldVector<double> p;
+// p = (*periodicMap)(coords[*it]);
+// TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
+//
+// association.push_back(std::make_pair(min(*it, idx2), max(*it, idx2)));
+// }
+// }
/**
* get al list of pairs that describes wich DegreeOfFreedom are assiciated in the given mesh
@@ -274,28 +216,51 @@ namespace details {
* The 'meshIndicator' describes one boundary of the mesh and to each DegreeOfFreedom
* the 'periodicMap' associates a second DegreeOfFreedom.
**/
+// inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
+// AbstractFunction<bool, WorldVector<double> >* meshIndicator,
+// AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
+// std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
+// {
+// std::set<DegreeOfFreedom> indices;
+// details::getBoundaryIndices(feSpace, meshIndicator, indices);
+// MSG_DBG("Nr of indices: %d\n", indices.size());
+//
+// association.clear();
+//
+// DOFVector<WorldVector<double> > coords(feSpace, "coords");
+// feSpace->getMesh()->getDofIndexCoords(coords);
+// std::set<DegreeOfFreedom>::iterator it;
+// for (it = indices.begin(); it != indices.end(); it++)
+// {
+// DegreeOfFreedom idx2;
+// WorldVector<double> p;
+// p = (*periodicMap)(coords[*it]);
+// TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
+//
+// association.push_back(std::make_pair(min(*it, idx2), max(*it, idx2));
+// }
+//
+// MSG_DBG("Nr of associations: %d\n", association.size());
+// }
+
inline void getPeriodicAssociation(const FiniteElemSpace* feSpace,
- AbstractFunction<bool, WorldVector<double> >* meshIndicator,
- AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
- std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
+ std::set<DegreeOfFreedom> indices,
+ AbstractFunction<WorldVector<double>, WorldVector<double> >* periodicMap,
+ std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &association)
{
- std::set<DegreeOfFreedom> indices;
- details::getBoundaryIndices(feSpace, meshIndicator, indices);
- MSG_DBG("Nr of indices: %d\n", indices.size());
-
association.clear();
DOFVector<WorldVector<double> > coords(feSpace, "coords");
feSpace->getMesh()->getDofIndexCoords(coords);
std::set<DegreeOfFreedom>::iterator it;
+
+ DegreeOfFreedom idx_;
+ WorldVector<double> p;
for (it = indices.begin(); it != indices.end(); it++)
{
- DegreeOfFreedom idx2;
- WorldVector<double> p;
p = (*periodicMap)(coords[*it]);
- TEST_EXIT(coords.getDofIdxAtPoint(p,idx2))("periodic association not found!\n");
-
- association.push_back(std::make_pair(*it, idx2));
+ TEST_EXIT(coords.getDofIdxAtPoint(p,idx_))("periodic association not found!\n");
+ association.push_back(std::make_pair(std::min(*it, idx_), std::max(*it, idx_)));
}
MSG_DBG("Nr of associations: %d\n", association.size());
@@ -547,11 +512,14 @@ struct PeriodicBcData {
void addToList(const FiniteElemSpace *feSpace, std::vector<ManualPeriodicBC> &list)
{
- std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > associations;
+ std::set<DegreeOfFreedom> indices;
if (nr == 0)
- details::getPeriodicAssociation(feSpace, meshIndicator, periodicMap, associations);
+ details::getBoundaryIndices(feSpace, meshIndicator, indices);
else
- details::getPeriodicAssociation(feSpace, nr, meshIndicator, periodicMap, associations);
+ details::getBoundaryIndices(feSpace, nr, meshIndicator, indices);
+
+ std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > associations;
+ details::getPeriodicAssociation(feSpace, indices, periodicMap, associations);
list.push_back(ManualPeriodicBC(row, associations));
}
int row;
diff --git a/extensions/base_problems/NavierStokesPhase_TaylorHood.cc b/extensions/base_problems/NavierStokesPhase_TaylorHood.cc
index 302ce2065542502069d446fbcf56489fdbafe773..8cce839c8ef0591bd23f8cbddeb7ec9011d37cf5 100644
--- a/extensions/base_problems/NavierStokesPhase_TaylorHood.cc
+++ b/extensions/base_problems/NavierStokesPhase_TaylorHood.cc
@@ -10,7 +10,8 @@ NavierStokesPhase_TaylorHood::NavierStokesPhase_TaylorHood(const std::string &na
beta(1.0),
epsilon(1.e-1),
alpha(3.0),
- fileWriterPhase(NULL)
+ fileWriterPhase(NULL),
+ phaseOld(NULL)
{
for (int i = 0; i < bcDOF.getSize(); ++i)
bcDOF[i] = NULL;
@@ -24,8 +25,14 @@ NavierStokesPhase_TaylorHood::NavierStokesPhase_TaylorHood(const std::string &na
NavierStokesPhase_TaylorHood::~NavierStokesPhase_TaylorHood()
{
- if (fileWriterPhase != NULL)
+ if (fileWriterPhase != NULL) {
delete fileWriterPhase;
+ fileWriterPhase = NULL;
+ }
+ if (phaseOld != NULL) {
+ delete phaseOld;
+ phaseOld = NULL;
+ }
}
@@ -44,29 +51,13 @@ void NavierStokesPhase_TaylorHood::transferInitialSolution(AdaptInfo *adaptInfo)
fileWriterPhase = new FileWriter(name + "->phase->output", getFeSpace()->getMesh(), phase);
for (int i = 0; i < dow; i++)
- transformDOF(prob->getSolution()->getDOFVector(0), phase, prob->getSolution()->getDOFVector(0), new AMDiS::Mult<double>);
+ transformDOF(prob->getSolution()->getDOFVector(i), phase, prob->getSolution()->getDOFVector(i), new AMDiS::Mult<double>);
super::transferInitialSolution(adaptInfo);
phaseOld->interpol(phase);
}
-void NavierStokesPhase_TaylorHood::closeTimestep(AdaptInfo *adaptInfo)
-{
- super::closeTimestep(adaptInfo);
- phaseOld->interpol(phase);
- writeFiles(adaptInfo, false);
-}
-
-
-void NavierStokesPhase_TaylorHood::writeFiles(AdaptInfo *adaptInfo, bool force)
-{ FUNCNAME("NavierStokesPhase_TaylorHood::closeTimestep()");
-
- super::writeFiles(adaptInfo, force);
- fileWriterPhase->writeFiles(adaptInfo, false);
-}
-
-
void NavierStokesPhase_TaylorHood::fillOperators()
{ FUNCNAME("NavierStokesPhase_TaylorHood::fillOperators()");
@@ -82,40 +73,34 @@ void NavierStokesPhase_TaylorHood::fillOperators()
Operator *opTime = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
opTime->addTerm(new Phase_ZOT(phase, density));
prob->addMatrixOperator(*opTime, i, i, getInvTau());
+
/// < (1/tau)*u_i^old , psi >
Operator *opTimeOld = new Operator(prob->getFeSpace(i), prob->getFeSpace(i));
- opTimeOld->addTerm(new Phase_ZOT(phaseOld, density));
- opTimeOld->setUhOld(getOldSolution(i));
+ opTimeOld->addTerm(new FactorPhase_ZOT(phaseOld, getOldSolution(i), density));
prob->addVectorOperator(*opTimeOld, i, getInvTau());
-
+
/// < u^old*grad(u_i^old) , psi >
- Operator *opUGradU0 = new Operator(prob->getFeSpace(i),prob->getFeSpace(i));
- opUGradU0->addTerm(new WorldVecPhase_FOT(phase,vel, -density), GRD_PHI);
+ Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i));
+ opUGradU0->addTerm(new WorldVecPhase_FOT(phase, vel, density), GRD_PHI);
opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
- if (nonLinTerm == 0) {
- prob->addVectorOperator(*opUGradU0, i);
- } else {
- prob->addVectorOperator(*opUGradU0, i, &theta1);
+ prob->addVectorOperator(*opUGradU0, i);
+
+ /// < u*grad(u_i^old) , psi >
+ for (unsigned k = 0; k < dow; ++k) {
+ Operator *opUGradU1 = new Operator(getFeSpace(i), getFeSpace(k));
+ opUGradU1->addTerm(new VecAndPartialDerivative_ZOT(
+ phase,
+ prob->getSolution()->getDOFVector(i), k, density));
+ prob->addMatrixOperator(*opUGradU1, i, k);
}
- if (nonLinTerm == 1) {
- /// < u'*grad(u_i^old) , psi >
- for (unsigned j = 0; j < dow; ++j) {
- Operator *opUGradU1 = new Operator(prob->getFeSpace(i),prob->getFeSpace(i));
- opUGradU1->addTerm(new VecAndPartialDerivative_ZOT(
- phase,
- prob->getSolution()->getDOFVector(i), j, density));
- prob->addMatrixOperator(*opUGradU1, i, j, &theta);
- }
- } else if (nonLinTerm == 2) {
- /// < u^old*grad(u'_i) , psi >
- for(unsigned j = 0; j < dow; ++j) {
- Operator *opUGradU2 = new Operator(prob->getFeSpace(i),prob->getFeSpace(i));
- opUGradU2->addTerm(new Vec2ProductPartial_FOT(
- phase,
- prob->getSolution()->getDOFVector(j), j, density), GRD_PHI);
- prob->addMatrixOperator(*opUGradU2, i, i, &theta);
- }
+ /// < u^old*grad(u_i) , psi >
+ for(unsigned k = 0; k < dow; ++k) {
+ Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
+ opUGradU2->addTerm(new Vec2ProductPartial_FOT(
+ phase,
+ vel[k], k, density), GRD_PHI);
+ prob->addMatrixOperator(*opUGradU2, i, i);
}
/// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity
@@ -147,6 +132,10 @@ void NavierStokesPhase_TaylorHood::fillOperators()
opDivU2->addTerm(new PartialDerivative_ZOT(phase,i));
prob->addMatrixOperator(*opDivU2, dow, i);
}
+
+ Operator *opNull = new Operator(getFeSpace(dow), getFeSpace(dow));
+ opNull->addTerm(new Simple_ZOT(0.0));
+ prob->addMatrixOperator(*opNull, dow, dow);
}
@@ -226,3 +215,18 @@ void NavierStokesPhase_TaylorHood::fillBoundaryConditions()
TEST_EXIT(!bcCond)("implicit boundary condition specified twice, by bcFct and bcDOF. This is not allowed!\n");
}
}
+
+
+void NavierStokesPhase_TaylorHood::closeTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("NavierStokesPhase_TaylorHood::closeTimestep()");
+ phaseOld->interpol(phase);
+ super::closeTimestep(adaptInfo);
+}
+
+
+void NavierStokesPhase_TaylorHood::writeFiles(AdaptInfo *adaptInfo, bool force)
+{ FUNCNAME("NavierStokesPhase_TaylorHood::writeFiles()");
+
+ super::writeFiles(adaptInfo, force);
+ fileWriterPhase->writeFiles(adaptInfo, false);
+}
diff --git a/extensions/base_problems/NavierStokes_TaylorHood.cc b/extensions/base_problems/NavierStokes_TaylorHood.cc
index 50251354916775e5d642d5e46e18a2173c667724..34bf3b71c27f8880d226ae90c2cdad6970db4cad 100644
--- a/extensions/base_problems/NavierStokes_TaylorHood.cc
+++ b/extensions/base_problems/NavierStokes_TaylorHood.cc
@@ -5,33 +5,17 @@ using namespace AMDiS;
NavierStokes_TaylorHood::NavierStokes_TaylorHood(const std::string &name_) :
super(name_),
- forceDBC(false),
- initialVelocityIsSet(false),
laplaceType(0),
nonLinTerm(2),
oldTimestep(0.0),
viscosity(1.0),
density(1.0),
- c(0.0),
- theta(0.5),
- theta1(0.5),
- minusTheta1(-0.5),
+ theta(0.5+1.e-2),
+ theta1(1-theta),
+ minusTheta1(-theta1),
velocity(NULL),
- fileWriter(NULL),
- initialX(NULL),
- initialY(NULL)
+ fileWriter(NULL)
{
- // force the homogeniouse dirichlet BC if combination of dirichlet and neumann BC
- // are set and AMDiS can not handle this combination automatically
- Initfile::get(name + "->force dirichlet bc", forceDBC);
- if (forceDBC) {
- int numDirichletPoints = 0;
- Initfile::get("number of dirichlet points", numDirichletPoints);
- dirichletPoints.resize(numDirichletPoints);
- for (int i = 0; i < numDirichletPoints; i++)
- Initfile::get("dirichlet point[" + Helpers::toString(i) + "]",dirichletPoints[i]);
- }
-
// parameters for navier-stokes
Initfile::get(name + "->viscosity", viscosity);
Initfile::get(name + "->density", density);
@@ -54,7 +38,7 @@ NavierStokes_TaylorHood::NavierStokes_TaylorHood(const std::string &name_) :
} else
dimension.set(1.0);
-
+ oldSolution.resize(dow);
for (size_t i = 0; i < dow; i++)
oldSolution[i] = NULL;
}
@@ -62,11 +46,6 @@ NavierStokes_TaylorHood::NavierStokes_TaylorHood(const std::string &name_) :
NavierStokes_TaylorHood::~NavierStokes_TaylorHood()
{ FUNCNAME("NavierStokes_TaylorHood::~NavierStokes_TaylorHood()");
-
- if (initialVelocityIsSet) {
- delete initialX;
- delete initialY;
- }
if (velocity != NULL) {
delete velocity;
@@ -101,6 +80,8 @@ void NavierStokes_TaylorHood::initData()
void NavierStokes_TaylorHood::solveInitialProblem(AdaptInfo *adaptInfo)
{ FUNCNAME("NavierStokes_TaylorHood::solveInitialProblem()");
+
+ double c = 1.0;
int initialVelocity = 0, flowDirection = 1;
Initfile::get(name + "->initial velocity", initialVelocity);
Initfile::get(name + "->initial velocity->value", c);
@@ -122,10 +103,6 @@ void NavierStokes_TaylorHood::solveInitialProblem(AdaptInfo *adaptInfo)
initialFcts[1-flowDirection] = new ConstantFct(0.0);
} else
throw(std::runtime_error("Unknown initial velocity."));
-
- initialX = initialFcts[0];
- initialY = initialFcts[1];
- initialVelocityIsSet = true;
MSG("solve initial problem...\n");
for (unsigned i = 0; i < dow; ++i) {
@@ -138,74 +115,64 @@ void NavierStokes_TaylorHood::solveInitialProblem(AdaptInfo *adaptInfo)
void NavierStokes_TaylorHood::transferInitialSolution(AdaptInfo *adaptInfo)
{ FUNCNAME("NavierStokes_TaylorHood::transferInitialSolution()");
- calcVelocity();
-
for (int i = 0; i < dow; i++)
prob->setExactSolution(prob->getSolution()->getDOFVector(i), i);
+ calcVelocity();
for (size_t i = 0; i < dow; i++)
oldSolution[i]->copy(*prob->getSolution()->getDOFVector(i));
- fileWriter->writeFiles(adaptInfo, false);
- writeFiles(adaptInfo, false);
-
- // initial parameters for detecting mesh changes
- oldMeshChangeIdx= getMesh()->getChangeIndex();
+ super::transferInitialSolution(adaptInfo);
}
void NavierStokes_TaylorHood::fillOperators()
{ FUNCNAME("NavierStokes_TaylorHood::fillOperators()");
+ WorldVector<DOFVector<double>* > vel;
+ for (unsigned k = 0; k < dow; k++)
+ vel[k] = prob->getSolution()->getDOFVector(k);
+
// fill operators for prob
for (unsigned i = 0; i < dow; ++i) {
/// < (1/tau)*u'_i , psi >
Operator *opTime = new Operator(getFeSpace(i), getFeSpace(i));
opTime->addTerm(new Simple_ZOT(density));
prob->addMatrixOperator(*opTime, i, i, getInvTau(), getInvTau());
+
/// < (1/tau)*u_i^old , psi >
- opTime->setUhOld(prob->getSolution()->getDOFVector(i));
- prob->addVectorOperator(*opTime, i, getInvTau(), getInvTau());
+ Operator *opTimeOld = new Operator(getFeSpace(i), getFeSpace(i));
+ opTimeOld->addTerm(new Phase_ZOT(getOldSolution(i), density));
+ prob->addVectorOperator(*opTimeOld, i, getInvTau(), getInvTau());
+
/// < u^old*grad(u_i^old) , psi >
Operator *opUGradU0 = new Operator(getFeSpace(i), getFeSpace(i));
- opUGradU0->addTerm(new WorldVector_FOT(velocity, -density), GRD_PHI);
+ opUGradU0->addTerm(new WorldVec_FOT(vel, density), GRD_PHI);
opUGradU0->setUhOld(prob->getSolution()->getDOFVector(i));
- if (nonLinTerm == 0) {
- prob->addVectorOperator(*opUGradU0, i);
- } else if (abs(theta1) > DBL_TOL) {
- prob->addVectorOperator(*opUGradU0, i, &theta1, &theta1);
+ prob->addVectorOperator(*opUGradU0, i);
+
+ /// < u*grad(u_i^old) , psi >
+ for (unsigned k = 0; k < dow; ++k) {
+ Operator *opUGradU1 = new Operator(getFeSpace(i), getFeSpace(k));
+ opUGradU1->addTerm(new PartialDerivative_ZOT(
+ prob->getSolution()->getDOFVector(i), k, density));
+ prob->addMatrixOperator(*opUGradU1, i, k);
}
- if (nonLinTerm == 1) {
- /// < u'*grad(u_i^old) , psi >
- for (unsigned j = 0; j < dow; ++j) {
- Operator *opUGradU1 = new Operator(getFeSpace(i),getFeSpace(i));
- opUGradU1->addTerm(new PartialDerivative_ZOT(
- prob->getSolution()->getDOFVector(i), j, density));
- prob->addMatrixOperator(*opUGradU1, i, j, &theta, &theta);
- }
- } else if (nonLinTerm == 2) {
- /// < u^old*grad(u'_i) , psi >
- for(unsigned j = 0; j < dow; ++j) {
- Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
- opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
- prob->getSolution()->getDOFVector(j), j, density), GRD_PHI);
- prob->addMatrixOperator(*opUGradU2, i, i, &theta, &theta);
- }
+ /// < u^old*grad(u_i) , psi >
+ for(unsigned k = 0; k < dow; ++k) {
+ Operator *opUGradU2 = new Operator(getFeSpace(i),getFeSpace(i));
+ opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
+ vel[k], k, density), GRD_PHI);
+ prob->addMatrixOperator(*opUGradU2, i, i);
}
-// for (int j = 0; j < dow; ++j) {
-// Operator *opNull = new Operator(getFeSpace(i), getFeSpace(j));
-// opNull->addTerm(new Simple_ZOT(0.0));
-// prob->addMatrixOperator(*opNull, i, j);
-// }
-
- /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity
+ /// Diffusion-Operator
addLaplaceTerm(i);
/// < p , d_i(psi) >
- Operator *opGradP = new Operator(getFeSpace(i),getFeSpace(dow));
+ Operator *opGradP = new Operator(getFeSpace(i), getFeSpace(dow));
opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
prob->addMatrixOperator(*opGradP, i, dow);
@@ -220,39 +187,43 @@ void NavierStokes_TaylorHood::fillOperators()
/// div(u) = 0
for (unsigned i = 0; i < dow; ++i) {
/// < d_i(u'_i) , psi >
- Operator *opDivU = new Operator(getFeSpace(dow),getFeSpace(i));
+ Operator *opDivU = new Operator(getFeSpace(dow), getFeSpace(i));
opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
prob->addMatrixOperator(*opDivU, dow, i);
}
+
+ Operator *opNull = new Operator(getFeSpace(dow), getFeSpace(dow));
+ opNull->addTerm(new Simple_ZOT(0.0));
+ prob->addMatrixOperator(*opNull, dow, dow);
}
void NavierStokes_TaylorHood::addLaplaceTerm(int i)
{ FUNCNAME("NavierStokes_TaylorHood::addLaplaceTerm()");
- /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
- if (laplaceType == 1) {
- for (unsigned j = 0; j < dow; ++j) {
- Operator *opLaplaceUi1 = new Operator(getFeSpace(i), getFeSpace(j));
- opLaplaceUi1->addTerm(new MatrixIJ_SOT(j, i, viscosity));
- prob->addMatrixOperator(*opLaplaceUi1, i, j, &theta, &theta);
-
- if (abs(minusTheta1) > DBL_TOL) {
- opLaplaceUi1->setUhOld(prob->getSolution()->getDOFVector(j));
- prob->addVectorOperator(*opLaplaceUi1, i, &minusTheta1, &minusTheta1);
- }
+ /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
+ if (laplaceType == 1) {
+ for (unsigned j = 0; j < dow; ++j) {
+ Operator *opLaplaceUi1 = new Operator(getFeSpace(i), getFeSpace(j));
+ opLaplaceUi1->addTerm(new MatrixIJ_SOT(j, i, viscosity));
+ prob->addMatrixOperator(*opLaplaceUi1, i, j, &theta, &theta);
+
+ if (abs(minusTheta1) > DBL_TOL) {
+ opLaplaceUi1->setUhOld(prob->getSolution()->getDOFVector(j));
+ prob->addVectorOperator(*opLaplaceUi1, i, &minusTheta1, &minusTheta1);
}
}
-
- /// < alpha*grad(u'_i) , grad(psi) >
- Operator *opLaplaceUi = new Operator(getFeSpace(i), getFeSpace(i));
- opLaplaceUi->addTerm(new Simple_SOT(viscosity));
- prob->addMatrixOperator(*opLaplaceUi, i, i, &theta, &theta);
-
- if (abs(minusTheta1) > DBL_TOL) {
- opLaplaceUi->setUhOld(prob->getSolution()->getDOFVector(i));
- prob->addVectorOperator(*opLaplaceUi, i, &minusTheta1, &minusTheta1);
- }
+ }
+
+ /// < alpha*grad(u'_i) , grad(psi) >
+ Operator *opLaplaceUi = new Operator(getFeSpace(i), getFeSpace(i));
+ opLaplaceUi->addTerm(new Simple_SOT(viscosity));
+ prob->addMatrixOperator(*opLaplaceUi, i, i, &theta, &theta);
+
+ if (abs(minusTheta1) > DBL_TOL) {
+ opLaplaceUi->setUhOld(prob->getSolution()->getDOFVector(i));
+ prob->addVectorOperator(*opLaplaceUi, i, &minusTheta1, &minusTheta1);
+ }
}
@@ -269,6 +240,8 @@ void NavierStokes_TaylorHood::closeTimestep(AdaptInfo *adaptInfo)
calcVelocity();
for (size_t i = 0; i < dow; i++)
oldSolution[i]->copy(*prob->getSolution()->getDOFVector(i));
+
+ super::closeTimestep(adaptInfo);
}
diff --git a/extensions/base_problems/NavierStokes_TaylorHood.h b/extensions/base_problems/NavierStokes_TaylorHood.h
index ebd93685de842d751243f5c5170e28d6e9aa5835..49c5f465e1f4bb511e67f7e094cc2cea80ed9b44 100644
--- a/extensions/base_problems/NavierStokes_TaylorHood.h
+++ b/extensions/base_problems/NavierStokes_TaylorHood.h
@@ -72,7 +72,6 @@ public: // methods
protected: // variables
- AbstractFunction<double, WorldVector<double> > *initialX, *initialY;
WorldVector<DOFVector<double>*> initialVel;
std::vector<WorldVector<double> > dirichletPoints;
@@ -89,7 +88,6 @@ protected: // variables
}
bool forceDBC;
- bool initialVelocityIsSet;
int laplaceType;
int nonLinTerm;
@@ -97,7 +95,7 @@ protected: // variables
double oldTimestep;
double viscosity;
double density;
- double c;
+
double theta;
double theta1;
double minusTheta1;
diff --git a/extensions/base_problems/chns.h b/extensions/base_problems/chns.h
index bc6cf319833f723335bbc4ca7fd88297b39b8a56..7a466dfe127cb2091eef160e89ca542055a7afb9 100644
--- a/extensions/base_problems/chns.h
+++ b/extensions/base_problems/chns.h
@@ -5,6 +5,8 @@
#include "AMDiS.h"
+using namespace AMDiS;
+
/** \brief
* Abstract function for Cahn-Hilliard mobility
*/
diff --git a/extensions/demo/CMakeLists.txt b/extensions/demo/CMakeLists.txt
index d6f1a931ee54376bbf758453905beda9a25b9e42..51d138c4dcab17d8e9ab68755e0ec8aaa8e2c40e 100644
--- a/extensions/demo/CMakeLists.txt
+++ b/extensions/demo/CMakeLists.txt
@@ -29,17 +29,32 @@ endif(AMDIS_FOUND)
add_executable("drivenCavity_twophase" ${drivenCavity_twophase})
target_link_libraries("drivenCavity_twophase" ${BASIS_LIBS})
- set(ddFsi src/diffuseDomainFsi.cc)
- add_executable("ddFsi" ${ddFsi})
- target_link_libraries("ddFsi" ${BASIS_LIBS})
+# set(ddFsi src/diffuseDomainFsi.cc)
+# add_executable("ddFsi" ${ddFsi})
+# target_link_libraries("ddFsi" ${BASIS_LIBS})
set(movingMesh src/movingMesh.cc)
add_executable("movingMesh" ${movingMesh})
target_link_libraries("movingMesh" ${BASIS_LIBS})
- set(rosenbrockTest src/rosenbrock_test.cc)
- add_executable("rosenbrockTest" ${rosenbrockTest})
- target_link_libraries("rosenbrockTest" ${BASIS_LIBS})
+# set(rosenbrockTest src/rosenbrock_test.cc)
+# add_executable("rosenbrockTest" ${rosenbrockTest})
+# target_link_libraries("rosenbrockTest" ${BASIS_LIBS})
+
+
+ set(navierStokesDd src/navierStokes_diffuseDomain.cc)
+ add_executable("navierStokesDd" ${navierStokesDd})
+ target_link_libraries("navierStokesDd" ${BASIS_LIBS})
+
+
+ set(navierStokesDd2 src/navierStokes_diffuseDomain2.cc)
+ add_executable("navierStokesDd2" ${navierStokesDd2})
+ target_link_libraries("navierStokesDd2" ${BASIS_LIBS})
+
+
+ set(navierStokesDd3 src/navierStokes_diffuseDomain3.cc)
+ add_executable("navierStokesDd3" ${navierStokesDd3})
+ target_link_libraries("navierStokesDd3" ${BASIS_LIBS})
#create the output dir
diff --git a/extensions/demo/init/reinit.inc.2d b/extensions/demo/init/reinit.inc.2d
index 4d7d822ea427b31fcc1ae530635259e800a8e4ca..7d5f518c246660802c1b5a1907576c36b7283751 100644
--- a/extensions/demo/init/reinit.inc.2d
+++ b/extensions/demo/init/reinit.inc.2d
@@ -1,5 +1,5 @@
-reinit->tolerance: 1.e-4
+reinit->tolerance: 1.e-2
reinit->maximal number of iteration steps: 100
reinit->Gauss-Seidel iteration: 1
reinit->infinity value: 1.e8
-reinit->boundary initialization: 3
+reinit->boundary initialization: 3
\ No newline at end of file