Skip to content
Snippets Groups Projects
Commit 99cb41c0 authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

complete definition of stiffness matrix A

parent 5d525c21
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 185 additions and 28 deletions
src/dg.cc
+ 185
− 28
View file @ 99cb41c0
......@@ -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];
}
}
}
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment