Hello everybody,
I'm using the deal.ii library for assembling a FE-system of equations of
the heat equation. I'm using a locally refined grid with *hanging nodes*.
The deal.ii library is only used for assembling since I want to use the
system of equations as optimization-constraints within a PDE-constrained
optimization problem. The drawback of this approach is that I cannot make
use of the deal-functions handling hanging nodes. This is due to the fact
that the hanging-node-constraints have to be considered within the system
of equations and cannot be applied afterwards with the
ConstraintMatrix::distribute(solution) function.
The upside is, the system of equation needn't to be positive definite,
since the NLP-solver I'm using handles every kind of equations. So, I
decided to integrate the hanging-node-constraints into the system of
equations by myself. From my point of view this is no big deal: I deleted
all rows considering hanging nodes of the system matrix and right hand side
and added the hanging-node-constraints instead.
For reasons of verification I solved this system of equations with a direct
solver and everything looked fine. But when I compared it to the solution
of the system of equations handling the hanging nodes with the
deal-functions, I recognized a difference between both solutions. So I went
on debugging by making the problem as small as possible: I changed the
dimension to 2D and generated a minimal grid:
9 --10 -- 2 ------ 3
| | | |
6 -- 7 -- 8 |
| | | |
4 -- 5 -- 0 ------ 1
having only one hanging node (8). Also the system of equations became very
easy:
(1 ) (1 )
( 1 ) (2 )
( ... ) x = (...)
( 1 ) (8 ) <-(constrained by node 0 and 2)
( 1 ) (9 )
( 1) (10 )
When I used my code:
void setup_system () {
dof_handler.distribute_dofs ( fe );
solution.reinit ( dof_handler.n_dofs() );
system_rhs.reinit ( dof_handler.n_dofs() );
constraints.clear ();
DoFTools::make_hanging_node_constraints ( dof_handler,
constraints );
constraints.close ();
DynamicSparsityPattern dsp ( dof_handler.n_dofs() );
DoFTools::make_sparsity_pattern ( dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ true );
sparsity_pattern.copy_from ( dsp );
system_matrix.reinit ( sparsity_pattern );
}
void assemble_system () {
FullMatrix<double> cell_matrix ( 1, 1 );
cell_matrix.set(0, 0, 1.0);
std::vector<types::global_dof_index> local_dof_indices ( 1 );
for (unsigned int i=0; i<dof_handler.n_dofs(); i++) {
local_dof_indices[0]= i;
user_distribute_local_to_global(cell_matrix, local_dof_indices,
system_matrix);
user_distribute_local_to_global(i, (double) i, system_rhs);
}
constrain_system_matrix();
}
void user_distribute_local_to_global(FullMatrix<double> cell_matrix,
std::vector<types::global_dof_index> local_dof_indices,
SparseMatrix<double> &global_matrix){
for ( unsigned int i = 0; i < cell_matrix.m(); i++ ) {
for ( unsigned int j = 0; j < cell_matrix.n(); j++ ) {
if ( !constraints.is_constrained ( local_dof_indices[i] ) ) {
global_matrix.set ( local_dof_indices[i], local_dof_indices[
j],
global_matrix.el ( local_dof_indices[i
], local_dof_indices[j] )
+ cell_matrix( i, j ) );
}
}
}
}
void user_distribute_local_to_global(unsigned int index, double value,
Vector<double> &global_vec){
if ( !constraints.is_constrained( index ) ) {
global_vec[index]=value;
}
}
void constrain_system_matrix() {
const std::vector< std::pair< types::global_dof_index, double > > *
constr_pair_ptr;
for ( unsigned int i=0; i<dof_handler.n_dofs(); i++ ) {
constr_pair_ptr=constraints.get_constraint_entries ( i );
if ( constr_pair_ptr ) {
for ( unsigned int j=0; j< ( *constr_pair_ptr ).size(); j++ ) {
system_matrix.set ( i, ( *constr_pair_ptr ) [j].first
,-1.0*( *constr_pair_ptr ) [j].second );
}
system_matrix.set(i,i, 1.);
}
}
}
- in my humble opinion - everything went fine. The solution is as expected:
x=(0, 1, 2, 3, 4, 5, 6, 7, 1, 9, 10)
^------------ (constrained by node 0 and 2: 0.5 *
0 + 0.5 * 2 = 1)
But when I used the deal-functions:
void setup_system () {
dof_handler.distribute_dofs ( fe );
solution.reinit ( dof_handler.n_dofs() );
system_rhs.reinit ( dof_handler.n_dofs() );
constraints.clear ();
DoFTools::make_hanging_node_constraints ( dof_handler,
constraints );
constraints.close ();
DynamicSparsityPattern dsp ( dof_handler.n_dofs() );
DoFTools::make_sparsity_pattern ( dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ true );
sparsity_pattern.copy_from ( dsp );
system_matrix.reinit ( sparsity_pattern );
}
void assemble_system () {
FullMatrix<double> cell_matrix ( 1, 1 );
cell_matrix.set(0, 0, 1.0);
std::vector<types::global_dof_index> local_dof_indices ( 1 );
for (unsigned int i=0; i<dof_handler.n_dofs(); i++) {
local_dof_indices[0]= i;
constraints.distribute_local_to_global ( cell_matrix,
local_dof_indices,
system_matrix);
constraints.distribute_local_to_global ( i,
(double) i,
system_rhs );
}
}
void solve () {
//...
constraints.distribute ( solution );
}
the solution confused me:
x=(2.33333, 1, 4.33333, 3, 4, 5, 6, 7, 3.33333, 9, 10)
Also the system of equations in the deal-functions-case looks wired:
(1.25 0.25 ) (4 )
( 1 ) (1 )
(0.25 1.25 ) (6 )
( 1 ) (3 )
( ... ) x = (...)
( 1 ) (0 )
( 1 ) (9 )
( 1) (10 )
So, what am I doing wrong? Am I using the deal-functions in a wrong manner?
In my humble opinion the way I'm treating the hanging node is correct. But
the fact that my results differ crucially from what deal is doing makes me
doubtful.
I'd be pleased about any help. Sorry for the long post.
Thanks a lot in advance,
Heinrich
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
