From 83d2b223af965b8dc6aa177baa39560152a7c465 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Fri, 17 Jan 2020 19:45:32 +0100
Subject: [PATCH] transfer assemblation of A and b to separate method
---
dune/mmdg/dg.hh | 740 ++++++++++++++++++++++++------------------------
1 file changed, 373 insertions(+), 367 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index e576d2d..e76112b 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -18,7 +18,7 @@ class DG
public:
using Scalar = typename GridView::ctype;
static constexpr int dim = GridView::dimension;
- using Matrix = Dune::DynamicMatrix<Scalar>;
+ using Matrix = Dune::DynamicMatrix<Scalar>; //NOTE: what is an appropriate sparse matrix type? -> BCRS
using Vector = Dune::DynamicVector<Scalar>;
//constructor
@@ -34,372 +34,8 @@ public:
const void operator() (const Scalar K, const Scalar mu)
{
- //NOTE: what is an appropriate sparse matrix type? -> BCRS
-
- //NOTE:
-// const int order = 3; //order of the quadrature rule
- for(const auto& elem : elements(gridView_))
- {
- const auto& geo = elem.geometry();
- std::cout << mapper_.index(elem) << "\t";
- for (int k = 0; k < geo.corners(); k++)
- {
- std::cout << geo.corner(k) << "\t";
- }
- std::cout << "\n\n";
- }
-
-
- //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 int elemIdx = mapper_.index(elem);
- const auto& geo = elem.geometry();
- const double elemVol = geo.volume();
- const auto& center = geo.center();
-
- std::cout << "===============================\nElement " << elemIdx;
- for (int k = 0; k < geo.corners(); k++)
- {
- std::cout << "\t" << geo.corner(k);
- }
- std::cout << "\n\n";
-
- //in the system of linear equations (SLE) Ad = b,
- //the index elemIdxSLE refers to the basis function phi_elem,0
- //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
- //basis function phi_elem,i
- const int elemIdxSLE = (dim + 1)*elemIdx;
-
- /* //TODO: can be done outside of the loop?
- const Dune::QuadratureRule<double,dim>& rule =
- Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
-
- //NOTE: how are quadrature rules in Dune applied correctly?
- for (const auto& ip : rule)
- {
- const auto qp = ip.position();
- //NOTE: is the volume of the element taken into account
- //automatically?
- const auto weight = ip.weight();
-
- //quadrature for int_elem q*phi_elem,0 dV
- b[elemIdxSLE] += weight * q(qp);
-
- //quadrature for int_elem q*phi_elem,i dV
- for (int i = 0; i < dim; i++)
- {
- b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
- }
- }
- */
-
- //NOTE: makeshift solution for source term q = -1
- b[elemIdxSLE] += -elemVol;
- for (int i = 0; i < dim; i++)
- {
- b[elemIdxSLE + i + 1] += -elemVol * center[i];
- }
-
- for (int i = 0; i < dim; i++)
- {
- //exact evaluation of
- // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
- A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
- }
-
- //iterate over all intersection with the boundary of elem
- for (const auto& intersct : intersections(gridView_, elem))
- {
- const auto& normal = intersct.centerUnitOuterNormal();
- const auto& intersctGeo = intersct.geometry();
- const double intersctVol = intersctGeo.volume();
- const auto& intersctCenter = intersctGeo.center();
-
- std::cout << "+++++++++++++++++++++++++++\nIntersection ";
- for (int k = 0; k < intersctGeo.corners(); k++)
- {
- std::cout << "\t" << intersctGeo.corner(k);
- }
- std::cout << "\tnormal: " << normal << "\n\n";
-
- //TODO: quadrature rule cannot be used for dim = 1!
-// const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
-// Dune::QuadratureRules<double,dim-1>::rule(
-// intersct.type(), 2);
-
- //storage for multiply used integral values,
- //note that only the lower diagonal and diagonal entries of
- //quadraticIntregrals are used
- Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
- Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
- //NOTE: is there a better type for symmetric matrix?
-
- for (int i = 0; i < dim; i++)
- {
- //use midpoint rule for exact evaluation of
- // int_intersct x_i ds
- linearIntegrals[i] = intersctVol * intersctCenter[i];
-
- for (int j = 0; j <= i; j++)
- {
- const auto& leftCorner = intersctGeo.corner(0);
- const auto& rightCorner = intersctGeo.corner(1);
- std::cout << "left corner: " << leftCorner << "\tright corner: " << rightCorner << "\n";
-
- quadraticIntregrals[i][j] = intersctVol / 3 *
- ( leftCorner[i] * leftCorner[j] +
- rightCorner[i] * rightCorner[j]
- + 0.5 * (leftCorner[i] * rightCorner[j] +
- leftCorner[j] * rightCorner[i]) );
-
- std::cout << i << ", " << j << "\t" << leftCorner[i] * leftCorner[j] << "\t" <<
- rightCorner[i] * rightCorner[j] << "\t" <<
- 0.5 * (leftCorner[i] * rightCorner[j] +
- leftCorner[j] * rightCorner[i]) << "\t" << intersctVol / 3 *
- ( leftCorner[i] * leftCorner[j] +
- rightCorner[i] * rightCorner[j]
- + 0.5 * (leftCorner[i] * rightCorner[j] +
- leftCorner[j] * rightCorner[i]) ) << "\t" << quadraticIntregrals[i][j] << "\n\n";
- /*
- //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;
- }
-*/
- //NOTE: probably unnecessary
- quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
- }
- }
- std::cout << "linearIntegrals:\n" << linearIntegrals << "\n\n";
- std::cout << "quadraticIntregrals:\n" << quadraticIntregrals << "\n\n";
-
- //exact evaluation of
- // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
- A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
- std::cout << elemIdxSLE << ", " << elemIdxSLE <<"\t" << intersctCenter<< ":\n"
- << mu * intersctVol << "\t" << A[elemIdx][elemIdx] << "\n\n";
-
- if (intersct.neighbor()) //intersct has neighboring element
- {
- std::cout << "------------------------------\nNeighbor\t";
- for (int k = 0; k < intersct.outside().geometry().corners(); k++)
- {
- std::cout << "\t" << intersct.outside().geometry().corner(k);
- }
- std::cout << "\n\n";
-
-
- //index of the neighboring element
- const int neighborIdx = mapper_.index(intersct.outside());
- const int neighborIdxSLE = (dim + 1)*neighborIdx;
-
- 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[elemIdxSLE + i + 1][elemIdxSLE] +=
- mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
- std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << ":\n"
- << mu * linearIntegrals[i] << "\t" <<
- - 0.5 * K * normal[i] * intersctVol << "\t"
- << A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
-
- 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[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
- += mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- + normal[j] * linearIntegrals[i]);
- std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j + 1 << ":\n"
- << mu * quadraticIntregrals[i][j] << "\t" <<
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- + normal[j] * linearIntegrals[i]) << "\t" <<
- A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
- }
- }
-
- 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[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
-
- //stiffness matrix A is symmetric
- A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];
-
- 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[elemIdxSLE + i + 1][neighborIdxSLE] +=
- -mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
- std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE << ":\n"
- << -mu * linearIntegrals[i] << "\t" <<
- 0.5 * K * normal[i] * intersctVol << "\t" <<
- A[elemIdxSLE + i + 1][neighborIdxSLE] << "\n\n";
-
- //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[elemIdxSLE][neighborIdxSLE + i + 1] +=
- -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
- std::cout << elemIdxSLE << ", " << neighborIdxSLE + i + 1<< ":\n"
- << -mu * linearIntegrals[i] << "\t" <<
- -0.5 * K * normal[i] * intersctVol << "\t" <<
- A[elemIdxSLE][neighborIdxSLE + i + 1] << "\n\n";
-
- //stiffness matrix A is symmetric
- A[neighborIdxSLE][elemIdxSLE + i + 1] +=
- A[elemIdxSLE + i + 1][neighborIdxSLE];
- A[neighborIdxSLE + i + 1][elemIdxSLE] +=
- A[elemIdxSLE][neighborIdxSLE + 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[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
- -mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[j] * linearIntegrals[i]
- - normal[i] * linearIntegrals[j]);
- std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE + j + 1<< ":\n"
- << -mu * quadraticIntregrals[i][j] << "\t" <<
- - 0.5 * K * (normal[j] * linearIntegrals[i]
- - normal[i] * linearIntegrals[j]) << "\t" <<
- A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] << "\n\n";
-
- //stiffness matrix A is symmetric
- A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
- A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
-
- if (i != 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[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
- -mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- - normal[j] * linearIntegrals[i]);
- std::cout << elemIdxSLE + j + 1 << ", " << neighborIdxSLE + i + 1<< ":\n"
- << -mu * quadraticIntregrals[i][j] << "\t" <<
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- - normal[j] * linearIntegrals[i]) << "\t" <<
- A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] << "\n\n";
-
- //stiffness matrix A is symmetric
- A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
- A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
- }
- }
- }
- }
- else //boundary facet
- {
- std::cout << "------------------------------\nBoundary\n\n";
- 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[elemIdxSLE + i + 1][elemIdxSLE] +=
- mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
- std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << "\t"
- << intersctCenter << ":\n" << mu * linearIntegrals[i] << "\t" <<
- - 0.5 * K * normal[i] * intersctVol << "\t"
- << A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
-
-
- 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[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
- mu * quadraticIntregrals[i][j]
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- + normal[j] * linearIntegrals[i]);
- std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j+ 1 << "\t"
- << intersctCenter << ":\n" << mu * quadraticIntregrals[i][j] << "\t" <<
- - 0.5 * K * (normal[i] * linearIntegrals[j]
- + normal[j] * linearIntegrals[i]) << "\t"
- << A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
- }
- }
- }
- }
-
- //stiffness matrix A is symmetric
- for (int i = 0; i < dim; i++)
- {
- A[elemIdxSLE][elemIdxSLE + i + 1] = A[elemIdxSLE + i + 1][elemIdxSLE];
-
- for(int j = 0; j < i; j++)
- {
- A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
- A[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
- }
- }
- }
-
- std::cout << A << std::endl;
+ //assemble stiffness matrix A and load vector b
+ assembleSLE(K, mu);
//NOTE: what would be an appropiate solver here?
A.solve(d, b);
@@ -488,6 +124,376 @@ public:
const int dof; //degrees of freedom
private:
+ //assemble stiffness matrix A and load vector b
+ void assembleSLE (const Scalar K, const Scalar mu)
+ {
+ //NOTE:
+//const int order = 3; //order of the quadrature rule
+
+ for (const auto& elem : elements(gridView_))
+ {
+ const auto& geo = elem.geometry();
+ std::cout << mapper_.index(elem) << "\t";
+ for (int k = 0; k < geo.corners(); k++)
+ {
+ std::cout << geo.corner(k) << "\t";
+ }
+ std::cout << "\n\n";
+ }
+
+ //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 int elemIdx = mapper_.index(elem);
+ const auto& geo = elem.geometry();
+ const double elemVol = geo.volume();
+ const auto& center = geo.center();
+
+ std::cout << "===============================\nElement " << elemIdx;
+ for (int k = 0; k < geo.corners(); k++)
+ {
+ std::cout << "\t" << geo.corner(k);
+ }
+ std::cout << "\n\n";
+
+ //in the system of linear equations (SLE) Ad = b,
+ //the index elemIdxSLE refers to the basis function phi_elem,0
+ //and the indices elemIdxSLE + i + 1, i = 1,...,dim, refer to the
+ //basis function phi_elem,i
+ const int elemIdxSLE = (dim + 1)*elemIdx;
+
+/* //TODO: can be done outside of the loop?
+ const Dune::QuadratureRule<double,dim>& rule =
+ Dune::QuadratureRules<double,dim>::rule(geo.type(),order);
+
+ //NOTE: how are quadrature rules in Dune applied correctly?
+ for (const auto& ip : rule)
+ {
+ const auto qp = ip.position();
+ //NOTE: is the volume of the element taken into account
+ //automatically?
+ const auto weight = ip.weight();
+
+ //quadrature for int_elem q*phi_elem,0 dV
+ b[elemIdxSLE] += weight * q(qp);
+
+ //quadrature for int_elem q*phi_elem,i dV
+ for (int i = 0; i < dim; i++)
+ {
+ b[elemIdxSLE + i + 1] += weight * qp[i] * q(qp);
+ }
+ }
+*/
+
+ //NOTE: makeshift solution for source term q = -1
+ b[elemIdxSLE] += -elemVol;
+ for (int i = 0; i < dim; i++)
+ {
+ b[elemIdxSLE + i + 1] += -elemVol * center[i];
+ }
+
+ for (int i = 0; i < dim; i++)
+ {
+ //exact evaluation of
+ // int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
+ A[elemIdxSLE + i + 1][elemIdxSLE + i + 1] += K * elemVol;
+ }
+
+ //iterate over all intersection with the boundary of elem
+ for (const auto& intersct : intersections(gridView_, elem))
+ {
+ const auto& normal = intersct.centerUnitOuterNormal();
+ const auto& intersctGeo = intersct.geometry();
+ const double intersctVol = intersctGeo.volume();
+ const auto& intersctCenter = intersctGeo.center();
+
+ std::cout << "+++++++++++++++++++++++++++\nIntersection ";
+ for (int k = 0; k < intersctGeo.corners(); k++)
+ {
+ std::cout << "\t" << intersctGeo.corner(k);
+ }
+ std::cout << "\tnormal: " << normal << "\n\n";
+
+ //TODO: quadrature rule cannot be used for dim = 1!
+// const Dune::QuadratureRule<double,dim-1>& secondOrderRule =
+// Dune::QuadratureRules<double,dim-1>::rule(
+// intersct.type(), 2);
+
+ //storage for multiply used integral values,
+ //note that only the lower diagonal and diagonal entries of
+ //quadraticIntregrals are used
+ Dune::FieldVector<Scalar, dim> linearIntegrals(0.0);
+ Dune::FieldMatrix<Scalar, dim, dim> quadraticIntregrals(0.0);
+ //NOTE: is there a better type for symmetric matrix?
+
+ for (int i = 0; i < dim; i++)
+ {
+ //use midpoint rule for exact evaluation of
+ // int_intersct x_i ds
+ linearIntegrals[i] = intersctVol * intersctCenter[i];
+
+ for (int j = 0; j <= i; j++)
+ {
+ const auto& leftCorner = intersctGeo.corner(0);
+ const auto& rightCorner = intersctGeo.corner(1);
+ std::cout << "left corner: " << leftCorner << "\tright corner: " << rightCorner << "\n";
+
+ quadraticIntregrals[i][j] = intersctVol / 3 *
+ ( leftCorner[i] * leftCorner[j] +
+ rightCorner[i] * rightCorner[j]
+ + 0.5 * (leftCorner[i] * rightCorner[j] +
+ leftCorner[j] * rightCorner[i]) );
+
+ std::cout << i << ", " << j << "\t" << leftCorner[i] * leftCorner[j] << "\t" <<
+ rightCorner[i] * rightCorner[j] << "\t" <<
+ 0.5 * (leftCorner[i] * rightCorner[j] +
+ leftCorner[j] * rightCorner[i]) << "\t" << intersctVol / 3 *
+ ( leftCorner[i] * leftCorner[j] +
+ rightCorner[i] * rightCorner[j]
+ + 0.5 * (leftCorner[i] * rightCorner[j] +
+ leftCorner[j] * rightCorner[i]) ) << "\t" << quadraticIntregrals[i][j] << "\n\n";
+/*
+ //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;
+ }
+*/
+ //NOTE: probably unnecessary
+ quadraticIntregrals[j][i] = quadraticIntregrals[i][j];
+ }
+ }
+ std::cout << "linearIntegrals:\n" << linearIntegrals << "\n\n";
+ std::cout << "quadraticIntregrals:\n" << quadraticIntregrals << "\n\n";
+
+ //exact evaluation of
+ // int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
+ A[elemIdxSLE][elemIdxSLE] += mu * intersctVol;
+ std::cout << elemIdxSLE << ", " << elemIdxSLE <<"\t" << intersctCenter<< ":\n"
+ << mu * intersctVol << "\t" << A[elemIdx][elemIdx] << "\n\n";
+
+ if (intersct.neighbor()) //intersct has neighboring element
+ {
+ std::cout << "------------------------------\nNeighbor\t";
+ for (int k = 0; k < intersct.outside().geometry().corners(); k++)
+ {
+ std::cout << "\t" << intersct.outside().geometry().corner(k);
+ }
+ std::cout << "\n\n";
+
+
+ //index of the neighboring element
+ const int neighborIdx = mapper_.index(intersct.outside());
+ const int neighborIdxSLE = (dim + 1)*neighborIdx;
+
+ 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[elemIdxSLE + i + 1][elemIdxSLE] +=
+ mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+ std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << ":\n"
+ << mu * linearIntegrals[i] << "\t" <<
+ - 0.5 * K * normal[i] * intersctVol << "\t"
+ << A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
+
+ 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[elemIdxSLE + i + 1][elemIdxSLE + j + 1]
+ += mu * quadraticIntregrals[i][j]
+ - 0.5 * K * (normal[i] * linearIntegrals[j]
+ + normal[j] * linearIntegrals[i]);
+ std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j + 1 << ":\n"
+ << mu * quadraticIntregrals[i][j] << "\t" <<
+ - 0.5 * K * (normal[i] * linearIntegrals[j]
+ + normal[j] * linearIntegrals[i]) << "\t" <<
+ A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
+ }
+ }
+
+ 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[elemIdxSLE][neighborIdxSLE] += -mu * intersctVol;
+
+ //stiffness matrix A is symmetric
+ A[neighborIdxSLE][elemIdxSLE] += A[elemIdxSLE][neighborIdxSLE];
+
+ 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[elemIdxSLE + i + 1][neighborIdxSLE] +=
+ -mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
+ std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE << ":\n"
+ << -mu * linearIntegrals[i] << "\t" <<
+ 0.5 * K * normal[i] * intersctVol << "\t" <<
+ A[elemIdxSLE + i + 1][neighborIdxSLE] << "\n\n";
+
+ //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[elemIdxSLE][neighborIdxSLE + i + 1] +=
+ -mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+ std::cout << elemIdxSLE << ", " << neighborIdxSLE + i + 1<< ":\n"
+ << -mu * linearIntegrals[i] << "\t" <<
+ -0.5 * K * normal[i] * intersctVol << "\t" <<
+ A[elemIdxSLE][neighborIdxSLE + i + 1] << "\n\n";
+
+ //stiffness matrix A is symmetric
+ A[neighborIdxSLE][elemIdxSLE + i + 1] +=
+ A[elemIdxSLE + i + 1][neighborIdxSLE];
+ A[neighborIdxSLE + i + 1][elemIdxSLE] +=
+ A[elemIdxSLE][neighborIdxSLE + 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[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
+ -mu * quadraticIntregrals[i][j]
+ - 0.5 * K * (normal[j] * linearIntegrals[i]
+ - normal[i] * linearIntegrals[j]);
+ std::cout << elemIdxSLE + i + 1 << ", " << neighborIdxSLE + j + 1<< ":\n"
+ << -mu * quadraticIntregrals[i][j] << "\t" <<
+ - 0.5 * K * (normal[j] * linearIntegrals[i]
+ - normal[i] * linearIntegrals[j]) << "\t" <<
+ A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1] << "\n\n";
+
+ //stiffness matrix A is symmetric
+ A[neighborIdxSLE + j + 1][elemIdxSLE + i + 1] +=
+ A[elemIdxSLE + i + 1][neighborIdxSLE + j + 1];
+
+ if (i != 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[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] +=
+ -mu * quadraticIntregrals[i][j]
+ - 0.5 * K * (normal[i] * linearIntegrals[j]
+ - normal[j] * linearIntegrals[i]);
+ std::cout << elemIdxSLE + j + 1 << ", " << neighborIdxSLE + i + 1<< ":\n"
+ << -mu * quadraticIntregrals[i][j] << "\t" <<
+ - 0.5 * K * (normal[i] * linearIntegrals[j]
+ - normal[j] * linearIntegrals[i]) << "\t" <<
+ A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1] << "\n\n";
+
+ //stiffness matrix A is symmetric
+ A[neighborIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+ A[elemIdxSLE + j + 1][neighborIdxSLE + i + 1];
+ }
+ }
+ }
+ }
+ else //boundary facet
+ {
+ std::cout << "------------------------------\nBoundary\n\n";
+ 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[elemIdxSLE + i + 1][elemIdxSLE] +=
+ mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
+ std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE << "\t"
+ << intersctCenter << ":\n" << mu * linearIntegrals[i] << "\t" <<
+ - 0.5 * K * normal[i] * intersctVol << "\t"
+ << A[elemIdxSLE + i + 1][elemIdxSLE] << "\n\n";
+
+
+ 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[elemIdxSLE + i + 1][elemIdxSLE + j + 1] +=
+ mu * quadraticIntregrals[i][j]
+ - 0.5 * K * (normal[i] * linearIntegrals[j]
+ + normal[j] * linearIntegrals[i]);
+ std::cout << elemIdxSLE + i + 1 << ", " << elemIdxSLE + j+ 1 << "\t"
+ << intersctCenter << ":\n" << mu * quadraticIntregrals[i][j] << "\t" <<
+ - 0.5 * K * (normal[i] * linearIntegrals[j]
+ + normal[j] * linearIntegrals[i]) << "\t"
+ << A[elemIdxSLE + i + 1][elemIdxSLE + j + 1] << "\n\n";
+ }
+ }
+ }
+ }
+
+ //stiffness matrix A is symmetric
+ for (int i = 0; i < dim; i++)
+ {
+ A[elemIdxSLE][elemIdxSLE + i + 1] = A[elemIdxSLE + i + 1][elemIdxSLE];
+
+ for(int j = 0; j < i; j++)
+ {
+ A[elemIdxSLE + j + 1][elemIdxSLE + i + 1] =
+ A[elemIdxSLE + i + 1][elemIdxSLE + j + 1];
+ }
+ }
+ }
+
+ std::cout << A << std::endl;
+ }
+
+
Matrix A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
--
GitLab