-
Notifications
You must be signed in to change notification settings - Fork 708
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add MatrixFreeTools::ElementBirthAndDeathMatrixFree #14121
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,7 @@ | ||
New: The new class MatrixFreeTools::ElementActivationAndDeactivationMatrixFree | ||
is a wrapper around MatrixFree designed to deal with | ||
DoFHandler objects involving cells without degrees of freedom, i.e., | ||
cells using FENothing as element type. The class helps to implement the | ||
"element birth and death technique". | ||
<br> | ||
(Peter Munch, Sebastian Proell, 2022/07/09) |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -157,6 +157,188 @@ namespace MatrixFreeTools | |
const unsigned int first_selected_component = 0); | ||
|
||
|
||
|
||
/** | ||
* A wrapper around MatrixFree to help users to deal with DoFHandler | ||
* objects involving cells without degrees of freedom, i.e., | ||
* cells using FENothing as element type. In the following we call such cells | ||
* deactivated. All other cells are activated. In contrast to MatrixFree, | ||
* this class skips deactivated cells and faces between activated and | ||
* deactivated cells are treated as boundary faces. | ||
Comment on lines
+165
to
+167
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think it would be important to say that this class considers only the first There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I have added a comment. To lift this restriction is the next step; see #14121 (comment). |
||
*/ | ||
template <int dim, | ||
typename Number, | ||
typename VectorizedArrayType = VectorizedArray<Number>> | ||
class ElementActivationAndDeactivationMatrixFree | ||
{ | ||
public: | ||
/** | ||
* Struct that helps to configure | ||
* ElementActivationAndDeactivationMatrixFree. | ||
*/ | ||
struct AdditionalData | ||
{ | ||
/** | ||
* Constructor. | ||
*/ | ||
AdditionalData(const unsigned int dof_index = 0) | ||
: dof_index(dof_index) | ||
{} | ||
|
||
/** | ||
* Index of the DoFHandler within MatrixFree to be used. | ||
*/ | ||
unsigned int dof_index; | ||
}; | ||
|
||
/** | ||
* Reinitialize class based on a given MatrixFree instance. | ||
* Particularly, the index of the valid FE is determined. | ||
* | ||
* @note At the moment, only DoFHandler objects are accepted | ||
* that are initialized with FECollection objects with at most | ||
* two finite elements (i.e., `FENothing` and another finite | ||
* element). | ||
*/ | ||
void | ||
reinit(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free, | ||
const AdditionalData &additional_data = AdditionalData()) | ||
{ | ||
this->matrix_free = &matrix_free; | ||
|
||
std::vector<unsigned int> valid_fe_indices; | ||
|
||
const auto &fe_collection = | ||
matrix_free.get_dof_handler(additional_data.dof_index) | ||
.get_fe_collection(); | ||
|
||
for (unsigned int i = 0; i < fe_collection.size(); ++i) | ||
if (fe_collection[i].n_dofs_per_cell() > 0) | ||
valid_fe_indices.push_back(i); | ||
|
||
// TODO: relax this so that arbitrary number of valid | ||
// FEs can be accepted | ||
AssertDimension(valid_fe_indices.size(), 1); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Can we allow more than one There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Yes, we can. Also, I guess what you need in your case is to not loop over faces neighboring Before making the modifications I would like to get to a status where we agree that the class is useful, we found a name, and maybe merged the first version so that we can improve it from there. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Yes, this should only be optional from my perspective. For me this class can definitely be useful, currently I skip FENothing cells (or cell batches) manually. |
||
|
||
fe_index_valid = *valid_fe_indices.begin(); | ||
} | ||
|
||
/** | ||
* Loop over all activated cells. | ||
* | ||
* For the meaning of the parameters see MatrixFree::cell_loop(). | ||
*/ | ||
template <typename VectorTypeOut, typename VectorTypeIn> | ||
void | ||
cell_loop(const std::function<void( | ||
const MatrixFree<dim, Number, VectorizedArrayType> &, | ||
VectorTypeOut &, | ||
const VectorTypeIn &, | ||
const std::pair<unsigned int, unsigned int>)> &cell_operation, | ||
VectorTypeOut & dst, | ||
const VectorTypeIn & src, | ||
const bool zero_dst_vector = false) | ||
{ | ||
const auto ebd_cell_operation = [&](const auto &matrix_free, | ||
auto & dst, | ||
const auto &src, | ||
const auto range) { | ||
const auto category = matrix_free.get_cell_range_category(range); | ||
|
||
if (category != fe_index_valid) | ||
return; | ||
|
||
cell_operation(matrix_free, dst, src, range); | ||
}; | ||
|
||
matrix_free->template cell_loop<VectorTypeOut, VectorTypeIn>( | ||
ebd_cell_operation, dst, src, zero_dst_vector); | ||
} | ||
|
||
/** | ||
* Loop over all activated cells and faces. Faces between activated and | ||
* deactivated cells are treated as boundary faces with the boundary ID | ||
* numbers::internal_face_boundary_id. | ||
* | ||
* For the meaning of the parameters see MatrixFree::cell_loop(). | ||
*/ | ||
template <typename VectorTypeOut, typename VectorTypeIn> | ||
void | ||
loop(const std::function< | ||
void(const MatrixFree<dim, Number, VectorizedArrayType> &, | ||
VectorTypeOut &, | ||
const VectorTypeIn &, | ||
const std::pair<unsigned int, unsigned int>)> &cell_operation, | ||
const std::function< | ||
void(const MatrixFree<dim, Number, VectorizedArrayType> &, | ||
VectorTypeOut &, | ||
const VectorTypeIn &, | ||
const std::pair<unsigned int, unsigned int>)> &face_operation, | ||
const std::function< | ||
void(const MatrixFree<dim, Number, VectorizedArrayType> &, | ||
VectorTypeOut &, | ||
const VectorTypeIn &, | ||
const std::pair<unsigned int, unsigned int>, | ||
const bool)> &boundary_operation, | ||
VectorTypeOut & dst, | ||
const VectorTypeIn & src, | ||
const bool zero_dst_vector = false) | ||
{ | ||
const auto ebd_cell_operation = [&](const auto &matrix_free, | ||
auto & dst, | ||
const auto &src, | ||
const auto range) { | ||
const auto category = matrix_free.get_cell_range_category(range); | ||
|
||
if (category != fe_index_valid) | ||
return; | ||
|
||
cell_operation(matrix_free, dst, src, range); | ||
}; | ||
|
||
const auto ebd_internal_or_boundary_face_operation = | ||
[&](const auto &matrix_free, | ||
auto & dst, | ||
const auto &src, | ||
const auto range) { | ||
const auto category = matrix_free.get_face_range_category(range); | ||
|
||
const unsigned int type = | ||
static_cast<unsigned int>(category.first == fe_index_valid) + | ||
static_cast<unsigned int>(category.second == fe_index_valid); | ||
|
||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Would There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. That is invalid face and nothing is done and the face is skipped. I have added a comment. |
||
if (type == 0) | ||
return; // deactivated face -> nothing to do | ||
|
||
if (type == 1) // boundary face | ||
boundary_operation( | ||
matrix_free, dst, src, range, category.first == fe_index_valid); | ||
else if (type == 2) // internal face | ||
face_operation(matrix_free, dst, src, range); | ||
}; | ||
|
||
matrix_free->template loop<VectorTypeOut, VectorTypeIn>( | ||
ebd_cell_operation, | ||
ebd_internal_or_boundary_face_operation, | ||
ebd_internal_or_boundary_face_operation, | ||
dst, | ||
src, | ||
zero_dst_vector); | ||
} | ||
|
||
private: | ||
/** | ||
* Reference to the underlying MatrixFree object. | ||
*/ | ||
SmartPointer<const MatrixFree<dim, Number, VectorizedArrayType>> | ||
matrix_free; | ||
|
||
/** | ||
* Index of the valid FE. Currently, ony a single one is supported. | ||
*/ | ||
unsigned int fe_index_valid; | ||
}; | ||
|
||
// implementations | ||
|
||
#ifndef DOXYGEN | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this a term used anywhere else on the literature? If not, then I would suggest to with something like activation/ deactivation instead.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You mean
ElementActivationAndDeactivationMatrixFree
?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, something like that sounds better to me. Would that be acceptable?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Renamed!