Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Add AffineConstraint::add_constraint(). #16103

Merged
merged 4 commits into from
Oct 6, 2023
Merged

Add AffineConstraint::add_constraint(). #16103

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
24d606d
Add AffineConstraint::add_constraint().
bangerth Oct 5, 2023
ddbb7b9
Add a changelog entry.
bangerth Oct 5, 2023
866ffdc
Convert a bunch of places to use AffineConstraint::add_constraint().
bangerth Oct 5, 2023
54fa322
Convert some of the tutorials.
bangerth Oct 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
7 changes: 7 additions & 0 deletions doc/news/changes/minor/20231002Bangerth_2
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
New: The new function AffineConstraints::add_constraint() adds an
entire constraint all at once, rather than splitting the task across
AffineConstraint::add_line(), multiple calls to
AffineConstraint::add_entry(), and
AffineConstraint::set_inhomogeneity().
<br>
(Wolfgang Bangerth, 2023/10/05)
7 changes: 4 additions & 3 deletions examples/step-11/step-11.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,16 +151,17 @@ namespace Step11

// Then generate a constraints object with just this one constraint. First
// clear all previous content (which might reside there from the previous
// computation on a once coarser grid), then add this one line
// computation on a once coarser grid), then add this one constraint,
// constraining the <code>first_boundary_dof</code> to the sum of other
// boundary DoFs each with weight -1. Finally, close the constraints
// object, i.e. do some internal bookkeeping on it for faster processing
// of what is to come later:
mean_value_constraints.clear();
mean_value_constraints.add_line(first_boundary_dof);
std::vector<std::pair<types::global_dof_index, double>> rhs;
for (const types::global_dof_index i : boundary_dofs)
if (i != first_boundary_dof)
mean_value_constraints.add_entry(first_boundary_dof, i, -1);
rhs.emplace_back(i, -1.);
mean_value_constraints.add_constraint(first_boundary_dof, rhs);
mean_value_constraints.close();

// Next task is to generate a sparsity pattern. This is indeed a tricky
Expand Down
3 changes: 1 addition & 2 deletions examples/step-41/step-41.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -479,8 +479,7 @@ namespace Step41
0)
{
active_set.add_index(dof_index);
constraints.add_line(dof_index);
constraints.set_inhomogeneity(dof_index, obstacle_value);
constraints.add_constraint(dof_index, {}, obstacle_value);

solution(dof_index) = obstacle_value;

Expand Down
15 changes: 5 additions & 10 deletions examples/step-46/doc/intro.dox
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -494,21 +494,16 @@ for (const auto &cell: dof_handler.active_cell_iterators())
cell->face(f)->get_dof_indices (local_face_dof_indices, 0);
for (unsigned int i=0; i<local_face_dof_indices.size(); ++i)
if (stokes_fe.face_system_to_component_index(i).first < dim)
constraints.add_line (local_face_dof_indices[i]);
constraints.add_constraint (local_face_dof_indices[i], {}, 0.);
}
}
@endcode

The call <code>constraints.add_line(t)</code> tells the
AffineConstraints to start a new constraint for degree of freedom
<code>t</code> of the form $x_t=\sum_{l=0}^{N-1} c_{tl} x_l +
b_t$. Typically, one would then proceed to set individual coefficients
$c_{tl}$ to nonzero values (using AffineConstraints::add_entry) or set
$b_t$ to something nonzero (using
AffineConstraints::set_inhomogeneity); doing nothing as above, funny as
it looks, simply leaves the constraint to be $x_t=0$, which is exactly
The last line calls
AffineConstraints::add_constraint() and adds the constraint
<i>x<sub>local_face_dof_indices[i]</sub>=0</i>, which is exactly
what we need in the current context. The call to
FiniteElement::face_system_to_component_index makes sure that we only set
FiniteElement::face_system_to_component_index() makes sure that we only set
boundary values to zero for velocity but not pressure components.

Note that there are cases where this may yield incorrect results:
Expand Down
4 changes: 3 additions & 1 deletion examples/step-46/step-46.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -385,7 +385,9 @@ namespace Step46
++i)
if (stokes_fe.face_system_to_component_index(i).first <
dim)
constraints.add_line(local_face_dof_indices[i]);
constraints.add_constraint(local_face_dof_indices[i],
{},
0.);
}
}
}
Expand Down
70 changes: 59 additions & 11 deletions include/deal.II/lac/affine_constraints.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,9 +246,10 @@ namespace internal
const size_type index_in_constraint) const;

/**
* Return whether there is one row with indirect contributions (i.e.,
* there has been at least one constraint with non-trivial
* ConstraintLine).
* Return whether this object stores a constraint in which one degree
* of freedom is constrained against another degree of freedom (as
* opposed to having only constraints of the form $x_i=42$ where
* no other degree of freedom appears on the right side).
*/
bool
have_indirect_rows() const;
Expand Down Expand Up @@ -428,14 +429,7 @@ namespace internal
} // namespace internal
#endif

// TODO[WB]: We should have a function of the kind
// AffineConstraints::add_constraint (const size_type constrained_dof,
// const std::vector<std::pair<size_type, number> > &entries,
// const number inhomogeneity = 0);
// rather than building up constraints piecemeal through add_line/add_entry
// etc. This would also eliminate the possibility of accidentally changing
// existing constraints into something pointless, see the discussion on the
// mailing list on "Tiny bug in interpolate_boundary_values" in Sept. 2010.


/**
* This class implements dealing with linear (possibly inhomogeneous)
Expand Down Expand Up @@ -760,6 +754,41 @@ class AffineConstraints : public Subscriptor
* @{
*/

/**
* Add a constraint to this object. This function adds a constraint
* of the form
* @f[
* x_i = \sum_{j=1}^n w_j x_{k_j} + b
* @f]
* where $i$ is the number of the degree of freedom to be constrained
* and is provided by the @p constrained_dofs argument. The weights
* $w_j$ and indices $k_j$ are provided as pairs in the
* @p dependencies argument, and the inhomogeneity $b$ is
* provided by the last argument.
*
* As an example, if you want to add the constraint
* @f[
* x_{42} = 0.5 x_{12} + 0.5 x_{36} + 27
* @f]
* you would call this function as follows:
* @code
* constraints.add_constraint (42, {{12, 0.5}, {36, 0.5}}, 27.0);
* @endcode
* On the other hand, if (as one often wants to) you need a constraint
* of the kind
* @f[
* x_{42} = 27
* @f]
* you would call this function as follows:
* @code
* constraints.add_constraint (42, {}, 27.0);
* @endcode
*/
void
add_constraint(const size_type constrained_dof,
const std::vector<std::pair<size_type, number>> &dependencies,
const number inhomogeneity = 0);

/**
* Add a new line to the matrix. If the line already exists, then the
* function simply returns without doing anything.
Expand Down Expand Up @@ -2203,6 +2232,25 @@ inline AffineConstraints<number>::AffineConstraints(



template <typename number>
inline void
AffineConstraints<number>::add_constraint(
const size_type constrained_dof,
const std::vector<std::pair<size_type, number>> &dependencies,
const number inhomogeneity)
{
Assert(is_constrained(constrained_dof) == false,
ExcMessage("You cannot add a constraint for a degree of freedom "
"that is already constrained."));

add_line(constrained_dof);
for (const auto &dep : dependencies)
add_entry(constrained_dof, dep.first, dep.second);
Comment on lines +2247 to +2248
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In light of #16121, I believe this is wrong and should be add_entries(...). I will make a quick test to see if it works, otherwise I have to postpone it because my time is constrained.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have a follow-up patch that avoids the repeated calls. Hold on with your timings.

set_inhomogeneity(constrained_dof, inhomogeneity);
}



template <typename number>
inline void
AffineConstraints<number>::add_line(const size_type line_n)
Expand Down
97 changes: 39 additions & 58 deletions include/deal.II/numerics/vector_tools_boundary.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -519,11 +519,9 @@ namespace VectorTools
{
if (constraints.can_store_line(boundary_value.first) &&
!constraints.is_constrained(boundary_value.first))
{
constraints.add_line(boundary_value.first);
constraints.set_inhomogeneity(boundary_value.first,
boundary_value.second);
}
constraints.add_constraint(boundary_value.first,
{},
boundary_value.second);
}
}

Expand Down Expand Up @@ -564,11 +562,9 @@ namespace VectorTools
{
if (constraints.can_store_line(boundary_value.first) &&
!constraints.is_constrained(boundary_value.first))
{
constraints.add_line(boundary_value.first);
constraints.set_inhomogeneity(boundary_value.first,
boundary_value.second);
}
constraints.add_constraint(boundary_value.first,
{},
boundary_value.second);
}
}

Expand Down Expand Up @@ -955,16 +951,13 @@ namespace VectorTools
std::map<types::global_dof_index, number> boundary_values;
project_boundary_values(
mapping, dof, boundary_functions, q, boundary_values, component_mapping);
typename std::map<types::global_dof_index, number>::const_iterator
boundary_value = boundary_values.begin();
for (; boundary_value != boundary_values.end(); ++boundary_value)

for (const auto &boundary_value : boundary_values)
{
if (!constraints.is_constrained(boundary_value->first))
{
constraints.add_line(boundary_value->first);
constraints.set_inhomogeneity(boundary_value->first,
boundary_value->second);
}
if (!constraints.is_constrained(boundary_value.first))
constraints.add_constraint(boundary_value.first,
{},
boundary_value.second);
}
}

Expand Down Expand Up @@ -1934,16 +1927,12 @@ namespace VectorTools
face_dof_indices[dof]) &&
!(constraints.is_constrained(
face_dof_indices[dof])))
{
constraints.add_line(
face_dof_indices[dof]);
if (std::abs(dof_values[dof]) > 1e-13)
{
constraints.set_inhomogeneity(
face_dof_indices[dof],
dof_values[dof]);
}
}
constraints.add_constraint(
face_dof_indices[dof],
{},
(std::abs(dof_values[dof]) > 1e-13 ?
dof_values[dof] :
0));
}
}
}
Expand Down Expand Up @@ -2078,17 +2067,12 @@ namespace VectorTools
face_dof_indices[dof]) &&
!(constraints.is_constrained(
face_dof_indices[dof])))
{
constraints.add_line(
face_dof_indices[dof]);

if (std::abs(dof_values[dof]) > 1e-13)
{
constraints.set_inhomogeneity(
face_dof_indices[dof],
dof_values[dof]);
}
}
constraints.add_constraint(
face_dof_indices[dof],
{},
(std::abs(dof_values[dof]) > 1e-13 ?
dof_values[dof] :
0));
}
}
}
Expand Down Expand Up @@ -2232,12 +2216,11 @@ namespace VectorTools
cell->face_orientation(face),
cell->face_flip(face),
cell->face_rotation(face)))[first_vector_component])
{
constraints.add_line(face_dof_indices[i]);

if (std::abs(dof_values(i)) > 1e-14)
constraints.set_inhomogeneity(face_dof_indices[i], dof_values(i));
}
constraints.add_constraint(face_dof_indices[i],
{},
(std::abs(dof_values[i]) > 1e-14 ?
dof_values[i] :
0));
}

// dummy implementation of above function for all other dimensions
Expand Down Expand Up @@ -2529,12 +2512,11 @@ namespace VectorTools
for (unsigned int dof = 0; dof < n_dofs; ++dof)
if ((projected_dofs[dof] != 0) &&
!(constraints.is_constrained(dof)))
{
constraints.add_line(dof);

if (std::abs(dof_values[dof]) > 1e-14)
constraints.set_inhomogeneity(dof, dof_values[dof]);
}
constraints.add_constraint(dof,
{},
(std::abs(dof_values[dof]) > 1e-14 ?
dof_values[dof] :
0));

break;
}
Expand Down Expand Up @@ -2688,12 +2670,11 @@ namespace VectorTools
for (unsigned int dof = 0; dof < n_dofs; ++dof)
if ((projected_dofs[dof] != 0) &&
!(constraints.is_constrained(dof)))
{
constraints.add_line(dof);

if (std::abs(dof_values[dof]) > 1e-14)
constraints.set_inhomogeneity(dof, dof_values[dof]);
}
constraints.add_constraint(dof,
{},
(std::abs(dof_values[dof]) > 1e-14 ?
dof_values[dof] :
0));

break;
}
Expand Down
32 changes: 16 additions & 16 deletions include/deal.II/numerics/vector_tools_constraints.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1063,13 +1063,13 @@ namespace VectorTools
constraints.can_store_line(
same_dof_range[0]->first.dof_indices[i]))
{
const types::global_dof_index line =
dof_indices.dof_indices[i];
constraints.add_line(line);
if (std::fabs(b_values[i]) >
std::numeric_limits<double>::epsilon())
constraints.set_inhomogeneity(line, b_values[i]);
// no add_entries here
constraints.add_constraint(
dof_indices.dof_indices[i],
{},
(std::fabs(b_values[i]) >
std::numeric_limits<double>::epsilon() ?
b_values[i] :
0));
}

break;
Expand Down Expand Up @@ -1459,10 +1459,7 @@ namespace VectorTools
{
const Vector<double> b_value = dof_vector_to_b_values[dofs];
for (unsigned int d = 0; d < dim; ++d)
{
constraints.add_line(dofs[d]);
constraints.set_inhomogeneity(dofs[d], b_value(d));
}
constraints.add_constraint(dofs[d], {}, b_value(d));
continue;
}

Expand Down Expand Up @@ -1502,12 +1499,15 @@ namespace VectorTools
const unsigned int new_index = dofs[d];
if (!constraints.is_constrained(new_index))
{
constraints.add_line(new_index);
if (std::abs(normal[d]) > 1e-13)
constraints.add_entry(new_index,
dofs[constrained_index],
-normal[d]);
constraints.set_inhomogeneity(new_index, boundary_value[d]);
constraints.add_constraint(new_index,
{{dofs[constrained_index],
-normal[d]}},
boundary_value[d]);
else
constraints.add_constraint(new_index,
{},
boundary_value[d]);
}
}
}
Expand Down
6 changes: 3 additions & 3 deletions source/dofs/dof_tools_constraints.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -571,9 +571,9 @@ namespace DoFTools
if (std::fabs(face_constraints(row, i)) >= 1e-14 * abs_sum)
entries.emplace_back(primary_dofs[i],
face_constraints(row, i));
constraints.add_line(dependent_dofs[row]);
constraints.add_entries(dependent_dofs[row], entries);
constraints.set_inhomogeneity(dependent_dofs[row], 0.);
constraints.add_constraint(dependent_dofs[row],
entries,
/* inhomogeneity= */ 0.);
}
}

Expand Down