diff --git a/src/dg.cc b/src/dg.cc
index 559b3447a3a9b3ea801c019164269cbe7b64ee45..23d2a9b00cfb03197f8069c021aa19079090e9b1 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -67,10 +67,17 @@ int main(int argc, char** argv)
const auto q = [](const Coordinate& x) {return x.two_norm();}; //source term
const int order = 3; //order of the quadrature rule
+
+ //we use the basis
+ // phi_elem,0 (x) = indicator(elem);
+ // phi_elem,i (x) = x[i]*indicator(elem);
+ //for elem in elements(gridView) and i = 1,...,dim
for (const auto& elem : elements(gridView))
{
- const auto& geo = elem.geometry();
const int elemIdx = mapper.index(elem);
+ const auto& geo = elem.geometry();
+ const double elemVol = geo.volume();
+
//TODO: can be done outside of the loop?
const Dune::QuadratureRule<double,dim>& rule =
Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
@@ -81,60 +88,210 @@ int main(int argc, char** argv)
//NOTE: is the volume of the element taken into account automatically?
const auto weight = ip.weight();
- //basis function phi_{elem,0}(x) = indicator(elem);
+ //quadrature for int_elem q*phi_elem,0 dV
b[(1 + dim)*elemIdx] += weight * q(qp);
+
+ //quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim; i++)
{
- //basis function phi_{elem,i}(x) = x[i]*indicator(elem);
b[(1 + dim)*elemIdx + i + 1] += weight * qp[i] * q(qp);
}
}
for (int i = 0; i < dim; i++)
{
- //basis function phi_{elem,i}(x) = x[i]*indicator(elem);
- A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + i + 1] =
- K * geo.volume(); // \int_elem K\nabla\phi_i\cdot\nabla\phi_i dV
+ //exact evaluation of
+ // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + i + 1] = K * elemVol;
}
- //NOTE: Is there a neat way to evaluate the jump and average operators?
- // (loop over facets?)
- //NOTE: jump / average operators disappear as basis function are supported
- //on a single element?
-
//iterate over all intersection with the boundary of elem
for (const auto& intersct : intersections(gridView, elem))
{
- const auto& outerNormal = intersct.centerUnitOuterNormal();
+ const auto& normal = intersct.centerUnitOuterNormal();
const auto& intersctGeo = intersct.geometry();
+ const double intersctVol = intersctGeo.volume();
+ const auto& intersctCenter = intersctGeo.center();
- const Dune::QuadratureRule<double,dim-1>& rule =
- Dune::QuadratureRules<double,dim-1>::rule(intersct.type(),order);
+ const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
+ Dune::QuadratureRules<double,dim-1>::rule(intersct.type(), 2);
- A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] += mu * intersctGeo.volume();
+ //storage for multiply used integral values,
+ //note that only the lower diagonal and diagonal entries of
+ //quadraticIntregrals are used
+ Dune::FieldVector<double, dim> linearIntegrals(0.0);
+ Dune::FieldMatrix<double, dim, dim> quadraticIntregrals(0.0);
- for (const auto& ip : rule)
+ for (int i = 0; i < dim; i++)
{
- const auto qp = ip.position();
- const auto weight = ip.weight();
-
- //TODO
+ //use midpoint rule for exact evaluation of
+ // int_intersct x_i ds
+ linearIntegrals[i] = mu * intersctVol * 0.5 * intersctCenter[i];
+
+ for (int j = 0; j <= i; j++)
+ {
+ //use second order quadrature rule for exact evaluation of
+ // int_intersct x_i*x_j ds
+ //NOTE: use Simpson's rule instead manually?
+ for (const auto& ip : secondOrderRule)
+ {
+ const auto& qp = ip.position();
+ quadraticIntregrals[i][j] += //NOTE: volume necessary?
+ ip.weight() * qp[i] * qp[j] * intersctVol;
+ }
+ }
}
-
+ //exact evaluation of
+ // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
+ A[(1 + dim)*elemIdx][(1 + dim)*elemIdx] += mu * intersctGeo.volume();
if (intersct.neighbor()) //intersct has neighboring element
{
//index of the neighboring element
- const int neighborIndex = mapper.index(intersct.outside());
-
- // A[][] = 0.5* ...; //note that interior facets are considered twice
-
- //TODO
+ const int neighborIdx = mapper.index(intersct.outside());
+
+ for (int i = 0; i < dim; i++)
+ { //we use the relations
+ // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
+ // = mu * int_intersct x_i ds
+ //and
+ // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
+ // = 0.5 * K * normal[i] * vol(intersct)
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] =
+ mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+
+ //stiffness matrix A is symmetric
+ A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] =
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
+
+ for (int j = 0; j <= i; j++)
+ {
+ //we use the relations
+ // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
+ // = mu * int_intersct x_i*x_j ds
+ //and
+ // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
+ // = 0.5 * K * normal[i] * int_intersct x_j ds
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1] =
+ mu * quadraticIntregrals[i][j] - 0.5 * K * (normal[i] *
+ linearIntegrals[j] + normal[j] * linearIntegrals[i]);
+
+ //stiffness matrix A is symmetric
+ A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] =
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1];
+ }
+ }
+
+ if (neighborIdx > elemIdx)
+ { //make sure that each facet is considered only once
+ continue;
+ }
+
+ //exact evaluation of
+ // int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
+ A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx] = -mu * intersctVol;
+
+ //stiffness matrix A is symmetric
+ A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx] =
+ A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx];
+
+ for (int i = 0; i < dim; i++)
+ {
+ //we use the relations
+ // int_intersct mu*jump(phi_elem,i)*jump(phi_neighbor,0) ds
+ // = -mu*int_intersct x_i ds
+ //and
+ // int_intersct avg(K*grad(phi_neighbor,i))*jump(phi_elem,0) ds
+ // = 0.5 * K * normal[i] * vol(intersct)
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx] =
+ -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+
+ //we use the relations
+ // int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,i) ds
+ // = -mu*int_intersct x_i ds
+ //and
+ // int_intersct avg(K*grad(phi_neighbor,i))*jump(phi_elem,0) ds
+ // = 0.5 * K * normal[i] * vol(intersct)
+ A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1] =
+ -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+
+ //stiffness matrix A is symmetric
+ A[(1 + dim)*neighborIdx][(1 + dim)*elemIdx + i + 1] =
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx];
+ A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx] =
+ A[(1 + dim)*elemIdx][(1 + dim)*neighborIdx + i + 1];
+
+ for (int j = 0; j <= i; j++)
+ {
+ //we use the relations
+ // int_intersct mu*jump(phi_elem,i)*jump(phi_neighbor,j) ds
+ // = -mu*int_intersct x_i*x_j ds
+ //and
+ // int_intersct avg(K*grad(phi_neighbor,j))*jump(phi_elem,i) ds
+ // = 0.5 * K * normal[j] * int_intersct x_i ds
+ //as well as
+ // int_intersct avg(K*grad(phi_elem,i))*jump(phi_neighbor,j) ds
+ // = -0.5 * K * normal[i] * int_intersct x_j ds
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx + j + 1] =
+ -mu * quadraticIntregrals[i][j]
+ - 0.5 * K * (normal[j] * linearIntegrals[i]
+ - normal[i] * linearIntegrals[j]);
+
+ // int_intersct mu*jump(phi_elem,j)*jump(phi_neighbor,i) ds
+ // = -mu*int_intersct x_i*x_j ds
+ //and
+ // int_intersct avg(K*grad(phi_neighbor,i))*jump(phi_elem,j) ds
+ // = 0.5 * K * normal[i] * int_intersct x_j ds
+ //as well as
+ // int_intersct avg(K*grad(phi_elem,j))*jump(phi_neighbor,i) ds
+ // = -0.5 * K * normal[j] * int_intersct x_i ds
+ A[(1 + dim)*elemIdx + j + 1][(1 + dim)*neighborIdx + i + 1] =
+ -mu * quadraticIntregrals[i][j]
+ - 0.5 * K * (normal[i] * linearIntegrals[j]
+ - normal[j] * linearIntegrals[i]);
+
+ //stiffness matrix A is symmetric
+ A[(1 + dim)*neighborIdx + j + 1][(1 + dim)*elemIdx + i + 1] =
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*neighborIdx + j + 1];
+ A[(1 + dim)*neighborIdx + i + 1][(1 + dim)*elemIdx + j + 1] =
+ A[(1 + dim)*elemIdx + j + 1][(1 + dim)*neighborIdx + i + 1];
+ }
+ }
}
- else //boundary element
+ else //boundary facet
{
- //TODO
+ for (int i = 0; i < dim; i++)
+ { //we use the relations
+ // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,i) ds
+ // = mu * int_intersct x_i ds
+ //and for boundary facets
+ // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
+ // = K * normal[i] * vol(intersct)
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx] =
+ mu * linearIntegrals[i] - K * normal[i] * intersctVol;
+
+ //stiffness matrix A is symmetric
+ A[(1 + dim)*elemIdx][(1 + dim)*elemIdx + i + 1] =
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx];
+
+ for (int j = 0; j <= i; j++)
+ {
+ //we use the relations
+ // int_intersct mu*jump(phi_elem,i)*jump(phi_elem,j) ds
+ // = mu * int_intersct x_i*x_j ds
+ //and for boundary facets
+ // int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,j) ds
+ // = 0.5 * K * normal[i] * int_intersct x_j ds
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1] =
+ mu * quadraticIntregrals[i][j] - K * (normal[i] *
+ linearIntegrals[j] + normal[j] * linearIntegrals[i]);
+
+ //stiffness matrix A is symmetric
+ A[(1 + dim)*elemIdx + j + 1][(1 + dim)*elemIdx + i + 1] =
+ A[(1 + dim)*elemIdx + i + 1][(1 + dim)*elemIdx + j + 1];
+ }
+ }
}
}
}