OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

update indices for new block matrices

Browse Source
This commit is contained in:
Arne Morten Kvarving 2017-03-29 12:12:35 +02:00
parent 586ca0a871
commit f6e08a4048
2 changed files with 8 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
Apps/Common/CompatibleOperators.C
View File

@ -36,9 +36,9 @@ void CompatibleOperators::Weak::Convection(std::vector<Matrix>& EM,
double scale,
WeakOperators::ConvectionForm form)
{
static const double vidx[3][3] = {{1, 6, 7},
{10, 2, 11},
{14, 15, 3}};
static const double vidx[3][3] = {{ 1, 7, 8},
{12, 2, 13},
{17, 18, 3}};
size_t nsd = fe.grad(1).cols();
for (size_t m=1; m <= nsd; ++m)
for (size_t i=1; i <= fe.basis(m).size(); ++i)
@ -61,7 +61,7 @@ void CompatibleOperators::Weak::Gradient(std::vector<Matrix>& EM,
for (size_t n = 1; n <= nsd; ++n)
for (size_t i=1; i <= fe.basis(n).size(); ++i)
for (size_t j=1; j <= fe.basis(nsd+1).size(); ++j)
EM[8+4*(n-1)](i,j) += -scale*fe.grad(n)(i,n)*fe.basis(nsd+1)(j)*fe.detJxW;
EM[10+5*(n-1)](i,j) += -scale*fe.grad(n)(i,n)*fe.basis(nsd+1)(j)*fe.detJxW;
}
@ -86,7 +86,7 @@ void CompatibleOperators::Weak::Laplacian(std::vector<Matrix>& EM,
for (size_t i = 1; i <= fe.basis(m).size(); ++i)
for (size_t n = m; n <= nsd; n++)
for (size_t j = 1; j <= fe.basis(n).size(); ++j) {
int idx = m == n ? m : (m == 1 ? 5+n-m : 10+n-m);
int idx = m == n ? m : (m == 1 ? 6+n-m : 13);
EM[idx](i,j) += scale* fe.grad(n)(j,m) * fe.grad(m)(i,n) * fe.detJxW;
}
}

6
Apps/Common/Test/TestCompatibleOperators.C
View File

@ -145,12 +145,12 @@ TEST(TestCompatibleOperators, Laplacian)
check_matrix_equal(EM[2], EM_2_ref);
// stress formulation
std::vector<Matrix> EM_stress(21);
std::vector<Matrix> EM_stress(36);
EM_stress[1].resize(6,6);
EM_stress[7].resize(6,6);
EM_stress[12].resize(6,6);
EM_stress[2].resize(6,6);
EM_stress[6].resize(6,6);
EM_stress[10].resize(6,6);
CompatibleOperators::Weak::Laplacian(EM_stress, fe, 1.0, true);
std::cout << EM_stress[1] << std::endl;