Commit 7b9f8d5a authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

several changes concerning quadrature degrees for operators

parent 47183212
......@@ -23,7 +23,7 @@ using namespace std;
using namespace AMDiS;
Phase_SOT::Phase_SOT(DOFVectorBase<double>* phaseDV_, double fac_)
: SecondOrderTerm(6),
: SecondOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
phaseDV(phaseDV_),
fac(fac_)
{
......@@ -91,7 +91,7 @@ double Phase_SOT::f(const int iq) const
/* ----------------------------------------------------------- */
Phase_FOT::Phase_FOT(DOFVectorBase<double>* phaseDV_, double fac_)
: FirstOrderTerm(6),
: FirstOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
phaseDV(phaseDV_),
fac(fac_)
{
......@@ -147,7 +147,7 @@ WorldVector<double> Phase_FOT::f(const int iq) const
/* ----------------------------------------------------------- */
Phase_ZOT::Phase_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
: ZeroOrderTerm(6),
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
phaseDV(phaseDV_),
fac(fac_)
{
......@@ -197,7 +197,7 @@ double Phase_ZOT::f(const int iq) const
/* ----------------------------------------------------------- */
PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
: ZeroOrderTerm(6),
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
phaseDV(phaseDV_),
fac(fac_),
component(0)
......@@ -210,7 +210,7 @@ PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_)
}
PhaseInverse_ZOT::PhaseInverse_ZOT(DOFVectorBase<double>* phaseDV_, double fac_, WorldVector<double> *facVec_, int component_)
: ZeroOrderTerm(6),
: ZeroOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
phaseDV(phaseDV_),
fac(fac_),
component(component_),
......@@ -262,419 +262,8 @@ double PhaseInverse_ZOT::f(const int iq) const
/* ----------------------------------------------------------- */
Pow3_ZOT::Pow3_ZOT(DOFVectorBase<double> *rhoDV_, double fac_)
: ZeroOrderTerm(3),
rhoDV(rhoDV_),
fac(fac_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void Pow3_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void Pow3_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
void Pow3_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += f(iq) * fac;
}
void Pow3_ZOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double opFactor)
{
double factor = opFactor * fac;
for (int iq = 0; iq < nPoints; iq++)
result[iq] += factor * f(iq) * uhAtQP[iq];
}
double Pow3_ZOT::f(const int iq) const
{
return rho[iq]*sqr(rho[iq]);
}
/* ----------------------------------------------------------- */
Pow3Phase_ZOT::Pow3Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV, fac_),
rhoDV(rhoDV_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void Pow3Phase_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(elInfo,subAssembler,quad);
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void Pow3Phase_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(largeElInfo,smallElInfo,subAssembler,quad);
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
double Pow3Phase_ZOT::f(const int iq) const
{
return rho[iq]*sqr(rho[iq]);
}
/* ----------------------------------------------------------- */
Pow2_ZOT::Pow2_ZOT(DOFVectorBase<double> *rhoDV_, double fac_)
: ZeroOrderTerm(2),
rhoDV(rhoDV_),
fac(fac_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void Pow2_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void Pow2_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
void Pow2_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += f(iq) * fac;
}
void Pow2_ZOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double opFactor)
{
double factor = opFactor * fac;
for (int iq = 0; iq < nPoints; iq++)
result[iq] += factor * f(iq) * uhAtQP[iq];
}
double Pow2_ZOT::f(const int iq) const
{
return sqr(rho[iq]);
}
/* ----------------------------------------------------------- */
Pow2Phase_ZOT::Pow2Phase_ZOT(DOFVectorBase<double>* phaseDV, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV, fac_),
rhoDV(rhoDV_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void Pow2Phase_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(elInfo,subAssembler,quad);
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void Pow2Phase_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(largeElInfo,smallElInfo,subAssembler,quad);
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
double Pow2Phase_ZOT::f(const int iq) const
{
return sqr(rho[iq]);
}
/* ----------------------------------------------------------- */
ConstrainedFrac_ZOT::ConstrainedFrac_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV_,fac_),
rhoDV(rhoDV_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void ConstrainedFrac_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(elInfo, subAssembler,quad);
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void ConstrainedFrac_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
double ConstrainedFrac_ZOT::f(const int iq) const
{
return 1.0/std::max(1.e-6, rho[iq]+0.9);
}
/* ----------------------------------------------------------- */
ConstrainedFracSqr_ZOT::ConstrainedFracSqr_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double factor_)
: Phase_ZOT(phaseDV_,factor_),
rhoDV(rhoDV_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void ConstrainedFracSqr_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(elInfo, subAssembler,quad);
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void ConstrainedFracSqr_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
double ConstrainedFracSqr_ZOT::f(const int iq) const
{
return -1.0/std::max(1.e-6, sqr(rho[iq]+0.9));
}
/* ----------------------------------------------------------- */
ConstrainedPowImpl_ZOT::ConstrainedPowImpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_)
: ZeroOrderTerm(2),
rhoDV(rhoDV_),
p(p_),
rhoMin(rhoMin_),
fac(fac_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void ConstrainedPowImpl_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void ConstrainedPowImpl_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
void ConstrainedPowImpl_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += f(iq)*fac;
}
void ConstrainedPowImpl_ZOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double opFactor)
{
double factor = opFactor*fac;
for (int iq = 0; iq < nPoints; iq++)
result[iq] += factor * f(iq) * uhAtQP[iq];
}
double ConstrainedPowImpl_ZOT::f(const int iq) const
{
return 2.0*p*Helpers::toInt(0.0<=rhoMin-rho[iq]);
}
/* ----------------------------------------------------------- */
ConstrainedPowExpl_ZOT::ConstrainedPowExpl_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_)
: ZeroOrderTerm(2),
rhoDV(rhoDV_),
p(p_),
rhoMin(rhoMin_),
fac(fac_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void ConstrainedPowExpl_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void ConstrainedPowExpl_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
void ConstrainedPowExpl_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += f(iq)*fac;
}
void ConstrainedPowExpl_ZOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double opFactor)
{
double factor = opFactor*fac;
for (int iq = 0; iq < nPoints; iq++)
result[iq] += factor * f(iq) * uhAtQP[iq];
}
double ConstrainedPowExpl_ZOT::f(const int iq) const
{
return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
}
/* ----------------------------------------------------------- */
ConstrainedPow_ZOT::ConstrainedPow_ZOT(DOFVectorBase<double> *rhoDV_, double p_, double rhoMin_, double fac_)
: ZeroOrderTerm(2),
rhoDV(rhoDV_),
p(p_),
rhoMin(rhoMin_),
fac(fac_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void ConstrainedPow_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void ConstrainedPow_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
void ConstrainedPow_ZOT::getC(const ElInfo *elInfo, int nPoints, ElementVector &C)
{
for (int iq = 0; iq < nPoints; iq++)
C[iq] += f(iq)*fac;
}
void ConstrainedPow_ZOT::eval(int nPoints,
const mtl::dense_vector<double>& uhAtQP,
const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
const mtl::dense_vector<WorldMatrix<double> >& D2UhAtQP,
mtl::dense_vector<double>& result,
double opFactor)
{
double factor = opFactor*fac;
for (int iq = 0; iq < nPoints; iq++)
result[iq] += factor * f(iq) * uhAtQP[iq];
}
double ConstrainedPow_ZOT::f(const int iq) const
{
return -2.0*p*std::max(0.0, rhoMin-rho[iq]);
}
/* ----------------------------------------------------------- */
FactorPhase_ZOT::FactorPhase_ZOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *rhoDV_, double fac_)
: Phase_ZOT(phaseDV_,fac_),
rhoDV(rhoDV_)
{
TEST_EXIT(rhoDV_)("No vector rho!\n");
auxFeSpaces.insert(rhoDV_->getFeSpace());
}
void FactorPhase_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(elInfo, subAssembler,quad);
getVectorAtQPs(rhoDV, elInfo, subAssembler, quad, rho);
}
void FactorPhase_ZOT::initElement(const ElInfo* largeElInfo,
const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
Phase_ZOT::initElement(largeElInfo, smallElInfo, subAssembler,quad);
getVectorAtQPs(rhoDV, smallElInfo, largeElInfo, subAssembler, quad, rho);
}
double FactorPhase_ZOT::f(const int iq) const
{
return rho[iq];
}
/* ----------------------------------------------------------- */
WorldVecAndVecFct_FOT::WorldVecAndVecFct_FOT(WorldVector<DOFVectorBase<double>*> vecs_, DOFVectorBase<double> *phaseDV_, double fac_)
: FirstOrderTerm(2),
: FirstOrderTerm(vecs_[0]->getFeSpace()->getBasisFcts()->getDegree() + phaseDV_->getFeSpace()->getBasisFcts()->getDegree()),
phaseDV(phaseDV_),
fac(fac_)
{
......@@ -750,7 +339,7 @@ void WorldVecAndVecFct_FOT::eval(int nPoints,
/* ----------------------------------------------------------- */
WorldVec_FOT::WorldVec_FOT(WorldVector<DOFVector<double>*> vecs_, double fac_)
: FirstOrderTerm(2),
: FirstOrderTerm(vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()),
fac(fac_)
{
numVecs=vecs_.size();
......@@ -821,7 +410,7 @@ void WorldVec_FOT::eval(int nPoints,
/* ----------------------------------------------------------- */
WorldVector_FOT::WorldVector_FOT(DOFVector<WorldVector<double> >* dv, double fac_)
: FirstOrderTerm(1),
: FirstOrderTerm(dv->getFeSpace()->getBasisFcts()->getDegree()),
vec(dv),
fac(fac_)
{
......@@ -872,7 +461,7 @@ void WorldVector_FOT::eval(int nPoints,
/* ----------------------------------------------------------- */
WorldVecPhase_FOT::WorldVecPhase_FOT(DOFVectorBase<double> *phaseDV_, WorldVector<DOFVector<double>*> vecs_,double fac_)
: FirstOrderTerm(2),
: FirstOrderTerm(phaseDV_->getFeSpace()->getBasisFcts()->getDegree() + vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()),
phaseDV(phaseDV_),
fac(fac_)
{
......@@ -949,7 +538,7 @@ void WorldVecPhase_FOT::eval(int nPoints,
/* ----------------------------------------------------------- */
VecAndWorldVec_FOT::VecAndWorldVec_FOT(DOFVector<double>* vec0_, WorldVector<DOFVector<double>*> vecs_, AbstractFunction<double, double>* fct_, double fac_)
: FirstOrderTerm(2),
: FirstOrderTerm(fct_->getDegree() + vecs_[0]->getFeSpace()->getBasisFcts()->getDegree()),
vDV(vec0_),
fct(fct_),
fac(fac_)
......@@ -1033,7 +622,7 @@ void VecAndWorldVec_FOT::eval(int nPoints,
// vec1*vec2*d/dxi(...)
VecAndPartialDerivative_FOT::VecAndPartialDerivative_FOT(DOFVectorBase<double> *dv1,
int component_, double fac_)
: FirstOrderTerm(3),
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree()),
vec1DV(dv1),
fac(fac_),
component(component_)
......@@ -1087,7 +676,7 @@ VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<doubl
int i_,
int j_,
double fac_)
: FirstOrderTerm(4),
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree() + dv2->getFeSpace()->getBasisFcts()->getDegree() - 1),
vec1DV(dv1),
vec2DV(dv2),
fct(NULL),
......@@ -1107,7 +696,7 @@ VecAndPartialDerivativeIJ_FOT::VecAndPartialDerivativeIJ_FOT(DOFVectorBase<doubl
int j_,
BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, d_j(v2))
double fac_)
: FirstOrderTerm(4),
: FirstOrderTerm(fct_->getDegree()),
vec1DV(dv1),
vec2DV(dv2),
fct(fct_),
......@@ -1178,7 +767,7 @@ Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
DOFVectorBase<double> *dv2,
int component_,
double fac_)
: FirstOrderTerm(3),
: FirstOrderTerm(dv1->getFeSpace()->getBasisFcts()->getDegree() + dv2->getFeSpace()->getBasisFcts()->getDegree()),
vec1DV(dv1),
vec2DV(dv2),
fct(NULL),
......@@ -1198,7 +787,7 @@ Vec2AndPartialDerivative_FOT::Vec2AndPartialDerivative_FOT(
int component_,
BinaryAbstractFunction<double, double, double>* fct_, // = f(v1, v2)
double fac_)
: FirstOrderTerm(3),
: FirstOrderTerm(fct_->getDegree()),
vec1DV(dv1),
vec2DV(dv2),
fct(fct_),
......@@ -1260,229 +849,10 @@ void Vec2AndPartialDerivative_FOT::eval(int nPoints,
/* ----------------------------------------------------------- */
// vec1*vec2*d/dxi(...)
PhaseAndViscosityPartialDerivativeBoundary_FOT::PhaseAndViscosityPartialDerivativeBoundary_FOT(DOFVectorBase<double> *phaseDV_, DOFVectorBase<double> *chDV_,
int component_, int normalComponent_, double viscosity1_, double viscosity2_)
: FirstOrderTerm(3),
phaseDV(phaseDV_),
chDV(chDV_),
PartialDerivative_FOT::PartialDerivative_FOT(int component_, double fac_)
: FirstOrderTerm(0),
component(component_),
normalComponent(normalComponent_),
v1(viscosity1_),
v2(viscosity2_)
{
auxFeSpaces.insert(phaseDV->getFeSpace());
auxFeSpaces.insert(chDV->getFeSpace());
setB(component);
}
void PhaseAndViscosityPartialDerivativeBoundary_FOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(phaseDV, elInfo, subAssembler, quad, phase);
getVectorAtQPs(chDV, elInfo, subAssembler, quad, ch);
WorldVector<double> normal;
int side= -1;
for (int i= 0; i < elInfo->getMesh()->getDim()+1; ++i) {
if(elInfo->getBoundary(i) > 2) {
side= i;
break;
}
}
// elInfo->getNormal(side, normal);
normalJ = 1.0; //normal[normalComponent];
}
void PhaseAndViscosityPartialDerivativeBoundary_FOT::initElement(const ElInfo* largeElInfo, const ElInfo* smallElInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
getVectorAtQPs(phaseDV, smallElInfo, largeElInfo, subAssembler, quad, phase);
getVectorAtQPs(chDV, smallElInfo, largeElInfo, subAssembler, quad, ch);
WorldVector<double> normal;
int side= -1;
for (int i= 0; i < largeElInfo->getMesh()->getDim()+1; ++i) {
if(largeElInfo->getBoundary(i) > 2) {
side= i;