Commit 9a8273c2 authored by Santiago Ospina's avatar Santiago Ospina
Browse files

[Flux Reconstruction] Calculate piola transformation with entity coordinates

* Fix 'order' keyword in raviart thomas volume local basis
parent 3d862731
......@@ -684,25 +684,25 @@ public:
// x-direction fastest, y-direction middle, z-direction slowest
std::size_t global_iter = 0.0;
for(std::size_t i2=0; i2<order+1; ++i2)
for(std::size_t i1=0; i1<order+1; ++i1)
for(std::size_t i0=0; i0<order; ++i0) // reduced
for(std::size_t i2=0; i2<k+1; ++i2)
for(std::size_t i1=0; i1<k+1; ++i1)
for(std::size_t i0=0; i0<k; ++i0) // reduced
{
out[global_iter] = {poly_km1.p(i0,in[0])*poly_k.p(i1,in[1])*poly_k.p(i2,in[2]), 0.0, 0.0};
global_iter++;
}
for(std::size_t i2=0; i2<order+1; ++i2)
for(std::size_t i1=0; i1<order; ++i1) // reduced
for(std::size_t i0=0; i0<order+1; ++i0)
for(std::size_t i2=0; i2<k+1; ++i2)
for(std::size_t i1=0; i1<k; ++i1) // reduced
for(std::size_t i0=0; i0<k+1; ++i0)
{
out[global_iter] = {0.0, poly_k.p(i0,in[0])*poly_km1.p(i1,in[1])*poly_k.p(i2,in[2]), 0.0};
global_iter++;
}
for(std::size_t i2=0; i2<order; ++i2) // reduced
for(std::size_t i1=0; i1<order+1; ++i1)
for(std::size_t i0=0; i0<order+1; ++i0)
for(std::size_t i2=0; i2<k; ++i2) // reduced
for(std::size_t i1=0; i1<k+1; ++i1)
for(std::size_t i0=0; i0<k+1; ++i0)
{
out[global_iter] = {0.0, 0.0, poly_k.p(i0,in[0])*poly_k.p(i1,in[1])*poly_km1.p(i2,in[2])};
global_iter++;
......@@ -821,8 +821,8 @@ public:
return k;
}
private:
Dune::QkStuff::GaussLobattoLagrangePolynomials<DF,RF,order> poly_k;
Dune::QkStuff::GaussLobattoLagrangePolynomials<DF,RF,order-1> poly_km1;
Dune::QkStuff::GaussLobattoLagrangePolynomials<DF,RF,k> poly_k;
Dune::QkStuff::GaussLobattoLagrangePolynomials<DF,RF,k-1> poly_km1;
};
#endif /* Actual code */
......
......@@ -43,6 +43,10 @@ class RaviartThomasVolumeLocalFiniteElementMap
RaviartThomasVolumeLocalFiniteElement<DF,RF,k,d,gt>,d>
{
public:
/*-----------------------------------------------------------------------*//**
* @brief Constructor of RaviartThomasVolumeLocalFiniteElementMap
*/
RaviartThomasVolumeLocalFiniteElementMap() {}
#ifdef DOXYGEN /* Documentation of the code */
......
......@@ -473,17 +473,19 @@ public:
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(p_local_n,gradphiv_n);
// Piola transformation
auto BTransposed = gtface.jacobianTransposed(it.position());
auto detB = BTransposed.determinant();
auto BTransposed_s = ig.inside().geometry().jacobianTransposed(p_local_s);
auto detB_s = BTransposed_s.determinant();
for (unsigned int i=0; i<lfsv_s.size(); i++) {
Vector y_s(gradphiv_s[i]);
y_s /= detB;
BTransposed.mtv(y_s,gradphiv_s[i]);
y_s /= detB_s;
BTransposed_s.mtv(y_s,gradphiv_s[i]);
}
auto BTransposed_n = ig.inside().geometry().jacobianTransposed(p_local_n);
auto detB_n = BTransposed_n.determinant();
for (unsigned int i=0; i<lfsv_n.size(); i++) {
Vector y_n(gradphiv_n[i]);
y_n /= detB;
BTransposed.mtv(y_n,gradphiv_n[i]);
y_n /= detB_n;
BTransposed_n.mtv(y_n,gradphiv_n[i]);
}
dg_sign = -1.;
......@@ -711,12 +713,12 @@ public:
// (we assume lfsv = RTk)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
// Piola transformation
auto BTransposed = gtface.jacobianTransposed(it.position());
auto detB = BTransposed.determinant();
auto BTransposed_s = ig.inside().geometry().jacobianTransposed(p_local_s);
auto detB_s = BTransposed_s.determinant();
for (unsigned int i=0; i<lfsv_s.size(); i++) {
Vector y_s(gradphiv_s[i]);
y_s /= detB;
BTransposed.mtv(y_s,gradphiv_s[i]);
y_s /= detB_s;
BTransposed_s.mtv(y_s,gradphiv_s[i]);
}
dg_sign = -1.;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment