On 6/26/19 3:59 AM, Alexander wrote:
> indeed, this should be possible to do it that way. Would you mind putting an
> example to the mentioned github issue? I might give it a thought in the
> background if I see more details.
OK, attached is some code that builds the constraints that correspond to this
rather esoteric paper (but also one I enjoyed working on the most because it's
just so weird!):
https://www.math.colostate.edu/~bangerth/publications/2016-nonconforming.pdf
The paper identifies 3 different ways how one could define hanging node
constraints for this element, and they are implemented in the function I
attach. The code is a couple of years old, but if something doesn't compile
any more, it should be easy to update.
Best
Wolfgang
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/54c1affc-8964-b910-a137-dc88088355d8%40colostate.edu.
For more options, visit https://groups.google.com/d/optout.
template <int dim>
Solver<dim>::LinearSystem::
LinearSystem (const DoFHandler<dim> &dof_handler)
{
hanging_node_constraints.clear ();
// Implementation of the following code only for the 2d DSSY case
Assert (dim == 2, ExcNotImplemented());
Assert (dynamic_cast<const DSSY::FE_DSSY*>(&dof_handler.get_fe()),
ExcNotImplemented());
std::vector<types::global_dof_index> face_dof_indices(1);
// Identify the following strategies outlined at the end of
// section 4 of the paper:
//
// 1) All averages are equal
// 2) Average on coarse face equals average of averages on child
// faces
// 3) Average jump is zero along child edges
//
// The following implements these strategies. select which
// strategy you want here:
const unsigned int strategy = 3;
std::map<types::global_dof_index,std::set<types::global_dof_index> > depends_on;
for (typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell)
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->at_boundary(f) == false)
// check of the neighbor is refined; if so, we need to
// treat the constraints on this interface
if (cell->face(f)->has_children())
{
switch (strategy)
{
case 1:
{
cell->face(f)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_mother_face = face_dof_indices[0];
cell->face(f)->child(0)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_child_face_0 = face_dof_indices[0];
cell->face(f)->child(1)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_child_face_1 = face_dof_indices[0];
// now enter two constraints that each of the
// child face dofs are equal to the parent face
// dof
hanging_node_constraints.add_line (dof_on_child_face_0);
hanging_node_constraints.add_entry (dof_on_child_face_0,
dof_on_mother_face,
1);
hanging_node_constraints.add_line (dof_on_child_face_1);
hanging_node_constraints.add_entry (dof_on_child_face_1,
dof_on_mother_face,
1);
depends_on[dof_on_child_face_0].insert (dof_on_mother_face);
depends_on[dof_on_child_face_1].insert (dof_on_mother_face);
break;
}
case 2:
{
cell->face(f)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_mother_face = face_dof_indices[0];
cell->face(f)->child(0)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_child_face_0 = face_dof_indices[0];
cell->face(f)->child(1)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_child_face_1 = face_dof_indices[0];
// now enter a constraint that the parent face dof
// equals one half of the sum of the child face dofs
hanging_node_constraints.add_line (dof_on_mother_face);
hanging_node_constraints.add_entry (dof_on_mother_face,
dof_on_child_face_0, 0.5);
hanging_node_constraints.add_entry (dof_on_mother_face,
dof_on_child_face_1, 0.5);
depends_on[dof_on_mother_face].insert (dof_on_child_face_0);
depends_on[dof_on_mother_face].insert (dof_on_child_face_1);
break;
}
case 3:
{
cell->face(f)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_mother_face = face_dof_indices[0];
cell->face(f)->child(0)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_child_face_0 = face_dof_indices[0];
cell->face(f)->child(1)->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_child_face_1 = face_dof_indices[0];
// for each face, define what the left and right neighbors are
static const unsigned int left_neighbor_faces[4] = { 2, 3, 0, 1 };
static const unsigned int right_neighbor_faces[4] = { 3, 2, 1, 0 };
cell->face(left_neighbor_faces[f])->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_left_neighbor = face_dof_indices[0];
cell->face(right_neighbor_faces[f])->get_dof_indices (face_dof_indices);
const types::global_dof_index dof_on_right_neighbor = face_dof_indices[0];
// now enter a constraint that child face dofs are
// linear combinations of the parent face dofs
hanging_node_constraints.add_line (dof_on_child_face_0);
hanging_node_constraints.add_line (dof_on_child_face_1);
switch (f)
{
case 0:
case 2:
{
hanging_node_constraints.add_entry (dof_on_child_face_0,
dof_on_mother_face,
1);
hanging_node_constraints.add_entry (dof_on_child_face_0,
dof_on_left_neighbor,
1./4);
hanging_node_constraints.add_entry (dof_on_child_face_0,
dof_on_right_neighbor,
-1./4);
hanging_node_constraints.add_entry (dof_on_child_face_1,
dof_on_mother_face,
1);
hanging_node_constraints.add_entry (dof_on_child_face_1,
dof_on_left_neighbor,
-1./4);
hanging_node_constraints.add_entry (dof_on_child_face_1,
dof_on_right_neighbor,
+1./4);
break;
}
case 1:
case 3:
{
hanging_node_constraints.add_entry (dof_on_child_face_0,
dof_on_mother_face,
1);
hanging_node_constraints.add_entry (dof_on_child_face_0,
dof_on_left_neighbor,
-1./4);
hanging_node_constraints.add_entry (dof_on_child_face_0,
dof_on_right_neighbor,
+1./4);
hanging_node_constraints.add_entry (dof_on_child_face_1,
dof_on_mother_face,
1);
hanging_node_constraints.add_entry (dof_on_child_face_1,
dof_on_left_neighbor,
+1./4);
hanging_node_constraints.add_entry (dof_on_child_face_1,
dof_on_right_neighbor,
-1./4);
break;
}
default:
Assert (false, ExcInternalError());
}
depends_on[dof_on_child_face_0].insert (dof_on_mother_face);
depends_on[dof_on_child_face_0].insert (dof_on_left_neighbor);
depends_on[dof_on_child_face_0].insert (dof_on_right_neighbor);
depends_on[dof_on_child_face_1].insert (dof_on_mother_face);
depends_on[dof_on_child_face_1].insert (dof_on_left_neighbor);
depends_on[dof_on_child_face_1].insert (dof_on_right_neighbor);
break;
}
default:
Assert (false, ExcNotImplemented());
}
}
std::map<unsigned int, unsigned int> chain_lengths;
for (std::map<types::global_dof_index,std::set<types::global_dof_index> >::const_iterator
p = depends_on.begin();
p != depends_on.end(); ++p)
{
unsigned int chain_length = 0;
std::set<types::global_dof_index> this_set = p->second;
while (this_set.empty() == false)
{
++chain_length;
// for each dependency, check whether it is constrained itself
std::set<types::global_dof_index> next_set;
for (std::set<types::global_dof_index>::const_iterator q = this_set.begin();
q != this_set.end(); ++q)
if (depends_on.find(*q) != depends_on.end())
next_set.insert (depends_on[*q].begin(),
depends_on[*q].end());
this_set = next_set;
}
++chain_lengths[chain_length];
}
std::cout << "N_constraints: " << hanging_node_constraints.n_constraints() << std::endl;
std::cout << "Chain lengths: ";
for (std::map<unsigned int, unsigned int>::const_iterator cl = chain_lengths.begin();
cl != chain_lengths.end(); ++cl)
std::cout << cl->first << ":" << cl->second << ", ";
std::cout << std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
hanging_node_constraints.close ();
hanging_node_constraints.condense (dsp);
sparsity_pattern.copy_from(dsp);
matrix.reinit (sparsity_pattern);
rhs.reinit (dof_handler.n_dofs());
}
