Skip to content

Commit 46835b7

Browse files
knarfnitramknarfnitram
committed
Refactor is_dbc_applied into own function
and remove from interface
1 parent 8b4acef commit 46835b7

12 files changed

+189
-162
lines changed

‎src/core/linalg/src/sparse/4C_linalg_blocksparsematrix.cpp‎

Lines changed: 0 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -185,23 +185,6 @@ void Core::LinAlg::BlockSparseMatrixBase::apply_dirichlet(
185185
}
186186
}
187187

188-
/*----------------------------------------------------------------------*
189-
*----------------------------------------------------------------------*/
190-
bool Core::LinAlg::BlockSparseMatrixBase::is_dbc_applied(const Core::LinAlg::Map& dbcmap,
191-
bool diagonalblock, const Core::LinAlg::SparseMatrix* trafo) const
192-
{
193-
for (int rblock = 0; rblock < rows(); ++rblock)
194-
{
195-
for (int cblock = 0; cblock < cols(); ++cblock)
196-
{
197-
if (not matrix(rblock, cblock)
198-
.is_dbc_applied(dbcmap, diagonalblock and (rblock == cblock), trafo))
199-
return false;
200-
}
201-
}
202-
return true;
203-
}
204-
205188

206189
/*----------------------------------------------------------------------*
207190
*----------------------------------------------------------------------*/

‎src/core/linalg/src/sparse/4C_linalg_blocksparsematrix.hpp‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -92,9 +92,6 @@ namespace Core::LinAlg
9292

9393
void apply_dirichlet(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true) override;
9494

95-
/// derived
96-
bool is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true,
97-
const Core::LinAlg::SparseMatrix* trafo = nullptr) const override;
9895

9996
//@}
10097

‎src/core/linalg/src/sparse/4C_linalg_sparsematrix.cpp‎

Lines changed: 0 additions & 111 deletions
Original file line numberDiff line numberDiff line change
@@ -1580,117 +1580,6 @@ void Core::LinAlg::SparseMatrix::add_other(Core::LinAlg::BlockSparseMatrixBase&
15801580
FOUR_C_THROW("BlockSparseMatrix and SparseMatrix cannot be added");
15811581
}
15821582

1583-
/*----------------------------------------------------------------------*
1584-
*----------------------------------------------------------------------*/
1585-
bool Core::LinAlg::SparseMatrix::is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock,
1586-
const Core::LinAlg::SparseMatrix* trafo) const
1587-
{
1588-
if (not filled()) FOUR_C_THROW("The matrix must be filled!");
1589-
1590-
const int numdbcrows = dbcmap.num_my_elements();
1591-
const int* dbcrows = dbcmap.my_global_elements();
1592-
1593-
std::vector<int> gIndices(max_num_entries(), 0);
1594-
std::vector<int> gtIndices((trafo ? trafo->max_num_entries() : 0), 0);
1595-
1596-
bool isdbc = true;
1597-
1598-
for (int i = 0; i < numdbcrows; ++i)
1599-
{
1600-
const int row = dbcrows[i];
1601-
const int sys_rlid = row_map().lid(row);
1602-
1603-
// this can happen for blocks of a BlockSparseMatrix
1604-
if (sys_rlid == -1) continue;
1605-
1606-
int NumEntries = 0;
1607-
double* Values = nullptr;
1608-
int* Indices = nullptr;
1609-
extract_my_row_view(sys_rlid, NumEntries, Values, Indices);
1610-
1611-
std::fill(gIndices.begin(), gIndices.end(), 0.0);
1612-
for (int c = 0; c < NumEntries; ++c) gIndices[c] = col_map().gid(Indices[c]);
1613-
1614-
// handle a diagonal block
1615-
if (diagonalblock)
1616-
{
1617-
if (NumEntries == 0) FOUR_C_THROW("Row {} is empty and part of a diagonal block!", row);
1618-
1619-
if (trafo)
1620-
{
1621-
if (not trafo->filled()) FOUR_C_THROW("The trafo matrix must be filled!");
1622-
1623-
int tNumEntries = 0;
1624-
double* tValues = nullptr;
1625-
int* tIndices = nullptr;
1626-
1627-
const int trafo_rlid = trafo->row_map().lid(row);
1628-
trafo->epetra_matrix().ExtractMyRowView(trafo_rlid, tNumEntries, tValues, tIndices);
1629-
1630-
// get the global indices corresponding to the extracted local indices
1631-
std::fill(gtIndices.begin(), gtIndices.end(), 0.0);
1632-
for (int c = 0; c < tNumEntries; ++c) gtIndices[c] = trafo->col_map().gid(tIndices[c]);
1633-
1634-
for (int j = 0; j < tNumEntries; ++j)
1635-
{
1636-
int k = -1;
1637-
while (++k < NumEntries)
1638-
if (Indices[k] == tIndices[j]) break;
1639-
1640-
if (k == NumEntries)
1641-
FOUR_C_THROW("Couldn't find column index {} in row {}.", tIndices[j], row);
1642-
1643-
if (std::abs(Values[k] - tValues[j]) > std::numeric_limits<double>::epsilon())
1644-
{
1645-
isdbc = false;
1646-
break;
1647-
}
1648-
}
1649-
}
1650-
// handle standard diagonal blocks
1651-
// --> 1.0 on the diagonal
1652-
// --> 0.0 on all off-diagonals
1653-
else
1654-
{
1655-
for (int j = 0; j < NumEntries; ++j)
1656-
{
1657-
if (gIndices[j] != row and Values[j] > std::numeric_limits<double>::epsilon())
1658-
{
1659-
isdbc = false;
1660-
break;
1661-
}
1662-
else if (gIndices[j] == row)
1663-
if (std::abs(1.0 - Values[j]) > std::numeric_limits<double>::epsilon())
1664-
{
1665-
isdbc = false;
1666-
break;
1667-
}
1668-
}
1669-
}
1670-
}
1671-
// we expect only zeros on the off-diagonal blocks
1672-
else
1673-
{
1674-
for (int j = 0; j < NumEntries; ++j)
1675-
{
1676-
if (Values[j] > std::numeric_limits<double>::epsilon())
1677-
{
1678-
isdbc = false;
1679-
break;
1680-
}
1681-
}
1682-
}
1683-
1684-
// stop as soon as the initial status changed once
1685-
if (not isdbc) break;
1686-
}
1687-
1688-
int lisdbc = static_cast<int>(isdbc);
1689-
int gisdbc = 0;
1690-
gisdbc = Core::Communication::min_all(lisdbc, Core::Communication::unpack_epetra_comm(Comm()));
1691-
1692-
return (gisdbc == 1);
1693-
}
16941583

16951584
/*----------------------------------------------------------------------*
16961585
*----------------------------------------------------------------------*/

‎src/core/linalg/src/sparse/4C_linalg_sparsematrix.hpp‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -411,8 +411,6 @@ namespace Core::LinAlg
411411
* rows.
412412
*
413413
* */
414-
bool is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true,
415-
const Core::LinAlg::SparseMatrix* trafo = nullptr) const override;
416414

417415
//@}
418416

‎src/core/linalg/src/sparse/4C_linalg_sparseoperator.hpp‎

Lines changed: 1 addition & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -196,17 +196,7 @@ namespace Core::LinAlg
196196
/// manner.
197197
virtual void apply_dirichlet(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true) = 0;
198198

199-
/** \brief Return TRUE if all Dirichlet boundary conditions have been applied
200-
* to this matrix
201-
*
202-
* \param (in) dbcmap: DBC map holding all dbc dofs
203-
* \param (in) diagonalblock: Is this matrix a diagonalblock of a blocksparsematrix?
204-
* If it is only one block/matrix, this boolean should be TRUE.
205-
* \param (in) trafo: pointer to an optional trafo matrix (see LocSys).
206-
*
207-
* */
208-
virtual bool is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true,
209-
const Core::LinAlg::SparseMatrix* trafo = nullptr) const = 0;
199+
210200

211201
/// Returns the Epetra_Map object associated with the (full) domain of this operator.
212202
virtual const Map& domain_map() const = 0;

‎src/core/linalg/src/sparse/4C_linalg_utils_sparse_algebra_assemble.cpp‎

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -279,4 +279,161 @@ std::shared_ptr<Core::LinAlg::MapExtractor> Core::LinAlg::convert_dirichlet_togg
279279
return std::make_shared<Core::LinAlg::MapExtractor>(fullmap, dbcmap, freemap);
280280
}
281281

282+
bool Core::LinAlg::is_dirichlet_boundary_condition_already_applied(
283+
const Core::LinAlg::SparseMatrix& A, const Core::LinAlg::Map& dbcmap, bool diagonalblock,
284+
const std::shared_ptr<const Core::LinAlg::SparseMatrix>& trafo)
285+
{
286+
if (not A.filled()) FOUR_C_THROW("The matrix must be filled!");
287+
288+
const int numdbcrows = dbcmap.num_my_elements();
289+
const int* dbcrows = dbcmap.my_global_elements();
290+
291+
std::vector<int> gIndices(A.max_num_entries(), 0);
292+
std::vector<int> gtIndices((trafo ? trafo->max_num_entries() : 0), 0);
293+
294+
bool isdbc = true;
295+
296+
for (int i = 0; i < numdbcrows; ++i)
297+
{
298+
const int row = dbcrows[i];
299+
const int sys_rlid = A.row_map().lid(row);
300+
301+
// this can happen for blocks of a BlockSparseMatrix
302+
if (sys_rlid == -1) continue;
303+
304+
int NumEntries = 0;
305+
double* Values = nullptr;
306+
int* Indices = nullptr;
307+
A.extract_my_row_view(sys_rlid, NumEntries, Values, Indices);
308+
309+
std::fill(gIndices.begin(), gIndices.end(), 0.0);
310+
for (int c = 0; c < NumEntries; ++c) gIndices[c] = A.col_map().gid(Indices[c]);
311+
312+
// handle a diagonal block
313+
if (diagonalblock)
314+
{
315+
if (NumEntries == 0) FOUR_C_THROW("Row {} is empty and part of a diagonal block!", row);
316+
317+
if (trafo)
318+
{
319+
if (not trafo->filled()) FOUR_C_THROW("The trafo matrix must be filled!");
320+
321+
int tNumEntries = 0;
322+
double* tValues = nullptr;
323+
int* tIndices = nullptr;
324+
325+
const int trafo_rlid = trafo->row_map().lid(row);
326+
trafo->epetra_matrix().ExtractMyRowView(trafo_rlid, tNumEntries, tValues, tIndices);
327+
328+
// get the global indices corresponding to the extracted local indices
329+
std::fill(gtIndices.begin(), gtIndices.end(), 0.0);
330+
for (int c = 0; c < tNumEntries; ++c) gtIndices[c] = trafo->col_map().gid(tIndices[c]);
331+
332+
for (int j = 0; j < tNumEntries; ++j)
333+
{
334+
int k = -1;
335+
while (++k < NumEntries)
336+
if (Indices[k] == tIndices[j]) break;
337+
338+
if (k == NumEntries)
339+
FOUR_C_THROW("Couldn't find column index {} in row {}.", tIndices[j], row);
340+
341+
if (std::abs(Values[k] - tValues[j]) > std::numeric_limits<double>::epsilon())
342+
{
343+
isdbc = false;
344+
break;
345+
}
346+
}
347+
}
348+
// handle standard diagonal blocks
349+
// --> 1.0 on the diagonal
350+
// --> 0.0 on all off-diagonals
351+
else
352+
{
353+
for (int j = 0; j < NumEntries; ++j)
354+
{
355+
if (gIndices[j] != row and Values[j] > std::numeric_limits<double>::epsilon())
356+
{
357+
isdbc = false;
358+
break;
359+
}
360+
else if (gIndices[j] == row)
361+
if (std::abs(1.0 - Values[j]) > std::numeric_limits<double>::epsilon())
362+
{
363+
isdbc = false;
364+
break;
365+
}
366+
}
367+
}
368+
}
369+
// we expect only zeros on the off-diagonal blocks
370+
else
371+
{
372+
for (int j = 0; j < NumEntries; ++j)
373+
{
374+
if (Values[j] > std::numeric_limits<double>::epsilon())
375+
{
376+
isdbc = false;
377+
break;
378+
}
379+
}
380+
}
381+
382+
// stop as soon as the initial status changed once
383+
if (not isdbc) break;
384+
}
385+
386+
int lisdbc = static_cast<int>(isdbc);
387+
int gisdbc = 0;
388+
gisdbc = Core::Communication::min_all(lisdbc, Core::Communication::unpack_epetra_comm(A.Comm()));
389+
390+
return (gisdbc == 1);
391+
}
392+
393+
394+
395+
/*----------------------------------------------------------------------*
396+
*----------------------------------------------------------------------*/
397+
bool Core::LinAlg::is_dirichlet_boundary_condition_already_applied(
398+
const Core::LinAlg::SparseOperator& A, const Core::LinAlg::Map& dbcmap, bool diagonalblock,
399+
const std::shared_ptr<const Core::LinAlg::SparseMatrix>& trafo)
400+
{
401+
// Check if provided Operator corresponds to single sparse matrix
402+
auto const* sparse_matrix = dynamic_cast<const Core::LinAlg::SparseMatrix*>(&A);
403+
if (sparse_matrix != nullptr)
404+
{
405+
return Core::LinAlg::is_dirichlet_boundary_condition_already_applied(
406+
*sparse_matrix, dbcmap, diagonalblock, trafo);
407+
}
408+
// Or Sparse Block Matrix
409+
auto const* block = dynamic_cast<const Core::LinAlg::BlockSparseMatrixBase*>(&A);
410+
if (block != nullptr)
411+
{
412+
for (int rblock = 0; rblock < block->rows(); ++rblock)
413+
{
414+
for (int cblock = 0; cblock < block->cols(); ++cblock)
415+
{
416+
const Core::LinAlg::SparseMatrix& submat = block->matrix(rblock, cblock);
417+
418+
if (not Core::LinAlg::is_dirichlet_boundary_condition_already_applied(
419+
submat, dbcmap, diagonalblock && (rblock == cblock), trafo))
420+
{
421+
return false;
422+
}
423+
}
424+
}
425+
return true;
426+
}
427+
else
428+
{
429+
FOUR_C_THROW(
430+
"Unsupported SparseOperator type in "
431+
"is_dirichlet_boundary_condition_already_applied().");
432+
}
433+
434+
435+
436+
return false;
437+
}
438+
282439
FOUR_C_NAMESPACE_CLOSE

‎src/core/linalg/src/sparse/4C_linalg_utils_sparse_algebra_assemble.hpp‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -265,6 +265,34 @@ namespace Core::LinAlg
265265
std::shared_ptr<Core::LinAlg::MapExtractor> convert_dirichlet_toggle_vector_to_maps(
266266
const Core::LinAlg::Vector<double>& dbctoggle);
267267

268+
/*! \brief Return TRUE if all Dirichlet boundary conditions have been applied
269+
to this sparse matrix
270+
271+
\param A (in) : Sparse Matrix A
272+
\param dbcmap (in) : DBC map holding all dbc dofs
273+
\param bool (in) diagonalblock: Is this matrix a diagonalblock of a blocksparsematrix?
274+
If it is only one block/matrix, this boolean should be TRUE.
275+
\param trafo (in) : shared pointer to an optional trafo matrix (see LocSys).
276+
*/
277+
bool is_dirichlet_boundary_condition_already_applied(const Core::LinAlg::SparseMatrix& A,
278+
const Core::LinAlg::Map& dbcmap, bool diagonalblock,
279+
const std::shared_ptr<const Core::LinAlg::SparseMatrix>& trafo);
280+
281+
282+
/*! \brief Return TRUE if all Dirichlet boundary conditions have been applied
283+
to this sparse operator
284+
285+
\param A (in) : SparseOperator A
286+
\param dbcmap (in) : DBC map holding all dbc dofs
287+
\param bool (in) diagonalblock: Is this matrix a diagonalblock of a blocksparsematrix?
288+
If it is only one block/matrix, this boolean should be TRUE.
289+
\param trafo (in) : shared pointer to an optional trafo matrix (see LocSys).
290+
291+
*/
292+
bool is_dirichlet_boundary_condition_already_applied(const Core::LinAlg::SparseOperator& A,
293+
const Core::LinAlg::Map& dbcmap, bool diagonalblock,
294+
const std::shared_ptr<const Core::LinAlg::SparseMatrix>& trafo);
295+
268296
} // namespace Core::LinAlg
269297

270298
FOUR_C_NAMESPACE_CLOSE

‎src/fsi/src/partitioned/nonlinear_solver/4C_fsi_nox_jacobian.cpp‎

Lines changed: 0 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -96,11 +96,6 @@ void NOX::FSI::FSIMatrixFree::apply_dirichlet(const Core::LinAlg::Map& dbcmap, b
9696
FOUR_C_THROW("Not implemented");
9797
}
9898

99-
bool NOX::FSI::FSIMatrixFree::is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock,
100-
const Core::LinAlg::SparseMatrix* trafo) const
101-
{
102-
FOUR_C_THROW("Not implemented");
103-
}
10499

105100
const Core::LinAlg::Map& NOX::FSI::FSIMatrixFree::domain_map() const
106101
{

‎src/fsi/src/partitioned/nonlinear_solver/4C_fsi_nox_jacobian.hpp‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -74,9 +74,6 @@ namespace NOX
7474

7575
void apply_dirichlet(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true) override;
7676

77-
bool is_dbc_applied(const Core::LinAlg::Map& dbcmap, bool diagonalblock = true,
78-
const Core::LinAlg::SparseMatrix* trafo = nullptr) const override;
79-
8077
const Core::LinAlg::Map& domain_map() const override;
8178

8279
void add(const Core::LinAlg::SparseOperator& A, const bool transposeA, const double scalarA,

0 commit comments

Comments
 (0)