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

changed: optimize assembleL2Matrices

Browse Source 
- assemble to a dense element matrix
- use outer_product for mass matrix calculation
- use an array of element rhs vectors to allow use of Vector::add
This commit is contained in:
Arne Morten Kvarving
2017-11-09 00:27:19 +01:00
parent e07697c317
commit bb60197ece
4 changed files with 69 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
src/ASM/ASMs2Drecovery.C
View File

@@ -229,6 +229,8 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
// --- Integration loop over all Gauss points in each direction ----------
Matrix eA(p1*p2, p1*p2);
Vectors eB(sField.rows(), Vector(p1*p1));
int ip = (i2*ng1*nel1 + i1)*ng2;
for (int j = 0; j < ng2; j++, ip += ng1*(nel1-1))
for (int i = 0; i < ng1; i++, ip++)
@@ -246,19 +248,22 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the linear system A*x=B
for (size_t ii = 0; ii < phi.size(); ii++)
{
int inod = MNPC[iel][ii]+1;
for (size_t jj = 0; jj < phi.size(); jj++)
{
int jnod = MNPC[iel][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
// Integrate the mass matrix
eA.outer_product(phi, phi, true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi,sField(r,ip+1)*dJw);
}
for (int i = 0; i < p1*p2; ++i) {
for (int j = 0; j < p1*p2; ++j)
A(MNPC[iel][i]+1, MNPC[iel][j]+1) += eA(i+1, j+1);
int jp = MNPC[iel][i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
return true;

31
src/ASM/ASMs3Drecovery.C
View File

@@ -236,6 +236,8 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
// --- Integration loop over all Gauss points in each direction --------
Matrix eA(p1*p2*p3, p1*p2*p3);
Vectors eB(sField.rows(), Vector(p1*p2*p3));
int ip = ((i3*ng2*nel2 + i2)*ng1*nel1 + i1)*ng3;
for (int k = 0; k < ng3; k++, ip += ng2*(nel2-1)*ng1*nel1)
for (int j = 0; j < ng2; j++, ip += ng1*(nel1-1))
@@ -254,20 +256,23 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the linear system A*x=B
for (size_t ii = 0; ii < phi.size(); ii++)
{
int inod = MNPC[iel][ii]+1;
for (size_t jj = 0; jj < phi.size(); jj++)
{
int jnod = MNPC[iel][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
// Integrate the mass matrix
eA.outer_product(phi, phi, true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi,sField(r,ip+1)*dJw);
}
}
for (int i = 0; i < p1*p2*p3; ++i) {
for (int j = 0; j < p1*p2*p3; ++j)
A(MNPC[iel][i]+1, MNPC[iel][j]+1) += eA(i+1, j+1);
int jp = MNPC[iel][i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
return true;
}

29
src/ASM/LR/ASMu2Drecovery.C
View File

@@ -151,6 +151,8 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
phi.resize((**el).nBasisFunctions());
// --- Integration loop over all Gauss points in each direction ----------
Matrix eA(MNPC[iel-1].size(), MNPC[iel-1].size());
Vectors eB(sField.rows(), Vector(MNPC[iel-1].size()));
int ip = 0;
for (int j = 0; j < ng2; j++)
for (int i = 0; i < ng1; i++, ip++)
@@ -174,19 +176,22 @@ bool ASMu2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the linear system A*x=B
for (size_t ii = 0; ii < phi.size(); ii++)
{
int inod = MNPC[iel-1][ii]+1;
for (size_t jj = 0; jj < phi.size(); jj++)
{
int jnod = MNPC[iel-1][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
// Integrate the mass matrix
eA.outer_product(phi, phi, true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi,sField(r,ip+1)*dJw);
}
for (size_t i = 0; i < MNPC[iel-1].size(); ++i) {
for (size_t j = 0; j < MNPC[iel-1].size(); ++j)
A(MNPC[iel-1][i]+1, MNPC[iel-1][j]+1) += eA(i+1, j+1);
int jp = MNPC[iel-1][i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
return true;

29
src/ASM/LR/ASMu3Drecovery.C
View File

@@ -159,6 +159,8 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
phi.resize((**el).nBasisFunctions());
// --- Integration loop over all Gauss points in each direction ----------
Matrix eA(MNPC[iel-1].size(), MNPC[iel-1].size());
Vectors eB(sField.rows(), Vector(MNPC[iel-1].size()));
int ip = 0;
for (int k = 0; k < ng3; k++)
for (int j = 0; j < ng2; j++)
@@ -183,19 +185,22 @@ bool ASMu3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
if (dJw == 0.0) continue; // skip singular points
}
// Integrate the linear system A*x=B
for (size_t ii = 0; ii < phi.size(); ii++)
{
int inod = MNPC[iel-1][ii]+1;
for (size_t jj = 0; jj < phi.size(); jj++)
{
int jnod = MNPC[iel-1][jj]+1;
A(inod,jnod) += phi[ii]*phi[jj]*dJw;
}
for (size_t r = 1; r <= sField.rows(); r++)
B(inod+(r-1)*nnod) += phi[ii]*sField(r,ip+1)*dJw;
}
// Integrate the mass matrix
eA.outer_product(phi, phi, true, dJw);
// Integrate the rhs vector B
for (size_t r = 1; r <= sField.rows(); r++)
eB[r-1].add(phi,sField(r,ip+1)*dJw);
}
for (size_t i = 0; i < MNPC[iel-1].size(); ++i) {
for (size_t j = 0; j < MNPC[iel-1].size(); ++j)
A(MNPC[iel-1][i]+1, MNPC[iel-1][j]+1) += eA(i+1, j+1);
int jp = MNPC[iel-1][i]+1;
for (size_t r = 0; r < sField.rows(); r++, jp += nnod)
B(jp) += eB[r](1+i);
}
}
return true;