From e1a3041b01288d48cc0bea1b430ebb0fe3d01f5a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Wed, 22 Jan 2020 19:04:00 +0100
Subject: [PATCH] use pointer for A in dg.hh
---
dune/mmdg/dg.hh | 65 ++++++++++++++++++++++---------------------------
1 file changed, 29 insertions(+), 36 deletions(-)
diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index d115dfa..3e69390 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -32,7 +32,7 @@ public:
dof((1 + dim) * gridView.size(0))
{
//initialize stiffness matrix A, load vector b and solution vector d
- A = Matrix(dof, dof, 4, 0.1, Matrix::implicit);
+ A = std::make_shared<Matrix>(dof, dof, 10, 1, Matrix::implicit); //TODO
b = Vector(dof);
d = Vector(dof);
}
@@ -43,7 +43,7 @@ public:
assembleSLE(K, mu);
Dune::InverseOperatorResult result;
- Dune::UMFPack<Matrix> solver(A);
+ Dune::UMFPack<Matrix> solver(*A);
solver.apply(d, b, result);
//storage for pressure data, for each element we store the
@@ -156,7 +156,7 @@ private:
{
//exact evaluation of
// int_elem K*grad(phi_elem,i)*grad(phi_elem,i) dV
- A.entry(elemIdxSLE + i + 1, elemIdxSLE + i + 1) += K * elemVol;
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + i + 1) += K * elemVol;
}
//iterate over all intersection with the boundary of elem
@@ -212,7 +212,7 @@ private:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_elem,0) ds
- A.entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
+ A->entry(elemIdxSLE, elemIdxSLE) += mu * intersctVol;
if (intersct.neighbor()) //intersct has neighboring element
{
@@ -227,7 +227,7 @@ private:
//and
// int_intersct avg(K*grad(phi_elem,i))*jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
- A.entry(elemIdxSLE + i + 1, elemIdxSLE) +=
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++)
@@ -238,7 +238,7 @@ private:
//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.entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1)
+= mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
@@ -252,11 +252,11 @@ private:
//exact evaluation of
// int_intersct mu*jump(phi_elem,0)*jump(phi_neighbor,0) ds
- A.entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
+ A->entry(elemIdxSLE, neighborIdxSLE) += -mu * intersctVol;
//stiffness matrix A is symmetric
- A.entry(neighborIdxSLE, elemIdxSLE) +=
- A.entry(elemIdxSLE, neighborIdxSLE);
+ A->entry(neighborIdxSLE, elemIdxSLE) +=
+ A->entry(elemIdxSLE, neighborIdxSLE);
for (int i = 0; i < dim; i++)
{
@@ -268,7 +268,7 @@ private:
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
- A.entry(elemIdxSLE + i + 1, neighborIdxSLE) +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE) +=
-mu * linearIntegrals[i] + 0.5 * K * normal[i] * intersctVol;
//we use the relations
@@ -279,14 +279,14 @@ private:
// int_intersct avg(K*grad(phi_neighbor,i))
// *jump(phi_elem,0) ds
// = 0.5 * K * normal[i] * vol(intersct)
- A.entry(elemIdxSLE, neighborIdxSLE + i + 1) +=
+ A->entry(elemIdxSLE, neighborIdxSLE + i + 1) +=
-mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
//stiffness matrix A is symmetric
- A.entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
- A.entry(elemIdxSLE + i + 1, neighborIdxSLE);
- A.entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
- A.entry(elemIdxSLE, neighborIdxSLE + i + 1);
+ A->entry(neighborIdxSLE, elemIdxSLE + i + 1) +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE);
+ A->entry(neighborIdxSLE + i + 1, elemIdxSLE) +=
+ A->entry(elemIdxSLE, neighborIdxSLE + i + 1);
for (int j = 0; j <= i; j++)
{
@@ -302,14 +302,14 @@ private:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_neighbor,j) ds
// = -0.5 * K * normal[i] * int_intersct x_j ds
- A.entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[j] * linearIntegrals[i]
- normal[i] * linearIntegrals[j]);
//stiffness matrix A is symmetric
- A.entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
- A.entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
+ A->entry(neighborIdxSLE + j + 1, elemIdxSLE + i + 1) +=
+ A->entry(elemIdxSLE + i + 1, neighborIdxSLE + j + 1);
if (i != j)
{
@@ -324,14 +324,14 @@ private:
// int_intersct avg(K*grad(phi_elem,j))
// *jump(phi_neighbor,i) ds
// = -0.5 * K * normal[j] * int_intersct x_i ds
- A.entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
+ A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
-mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
- normal[j] * linearIntegrals[i]);
//stiffness matrix A is symmetric
- A.entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
- A.entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
+ A->entry(neighborIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+ A->entry(elemIdxSLE + j + 1, neighborIdxSLE + i + 1);
}
}
}
@@ -346,7 +346,7 @@ private:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,0) ds
// = K * normal[i] * vol(intersct)
- A.entry(elemIdxSLE + i + 1, elemIdxSLE) +=
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE) +=
mu * linearIntegrals[i] - 0.5 * K * normal[i] * intersctVol;
for (int j = 0; j <= i; j++)
@@ -358,7 +358,7 @@ private:
// int_intersct avg(K*grad(phi_elem,i))
// *jump(phi_elem,j) ds
// = 0.5 * K * normal[i] * int_intersct x_j ds
- A.entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1) +=
mu * quadraticIntregrals[i][j]
- 0.5 * K * (normal[i] * linearIntegrals[j]
+ normal[j] * linearIntegrals[i]);
@@ -370,29 +370,22 @@ private:
//stiffness matrix A is symmetric
for (int i = 0; i < dim; i++)
{
- A.entry(elemIdxSLE, elemIdxSLE + i + 1) =
- A.entry(elemIdxSLE + i + 1, elemIdxSLE);
+ A->entry(elemIdxSLE, elemIdxSLE + i + 1) =
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE);
for(int j = 0; j < i; j++)
{
- A.entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
- A.entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
+ A->entry(elemIdxSLE + j + 1, elemIdxSLE + i + 1) =
+ A->entry(elemIdxSLE + i + 1, elemIdxSLE + j + 1);
}
}
}
- //NOTE: check if A is symmetric
- /* for (int i = 0; i < dof; i++)
- for (int j = 0; j < i; j++) */
- /*assert*//*if(std::abs(A[i][j] - A[j][i]) >
- std::numeric_limits<Scalar>::epsilon())
- std::cout << i << ", " << j << std::endl; */
-
- A.compress();
+ A->compress();
}
- Matrix A; //stiffness matrix
+ std::shared_ptr<Matrix> A; //stiffness matrix
Vector b; //load vector
Vector d; //solution vector
};
--
GitLab