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

use pointer for A in dg.hh

parent 295dac6f
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 29 additions and 36 deletions
dune/mmdg/dg.hh
+ 29
− 36
View file @ e1a3041b
......@@ -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
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment