Hi Andrew,
thanks for your reply. I also did think that the order might be a
problem, so I applied the periodic boundaries first as I wanted the
Dirichlet conditions to "win".
I now tried also your decision using two constraints and merging them
and set the dirichlet ones as the winner. But didn't hasn't had any
effect. I got the same matrix and rhs and the same result.
I than tried it also the other way round, which again gave the same
behavior.
So, something else must be the problem. In addition I found out the
number of dofs common to both constraints and explicitly skipped setting
periodic constraints to them. Still the same...
Btw, this behavior only occurs for one of the two components.
Timo
Andrew McBride wrote:
Hi Timo
In your implementation you only have one set of constraints. you apply the periodic conditions first and then the Dirichlet. The nodes common to the Dirichlet and periodic boundaries would then take the value of the Dirichlet condition (actually I don't think the order is the issue here, it's the fact that you constrain these nodes twice). What value do you want these nodes to take?
I do something similar in 3D but I merge different sets of constraints and
specify which set wins to get the final constraints.
does this help?
Andrew
On 14.04.2011, at 11:36, Timo Koellner wrote:
Hi folks,
I'm still dealing with the following problem (see "Coupled Set of PDE" March,
30th)
\nabla n = 0
\nabla \vec{E} = -e(n - \Theta(x))
In this first order form, there is the problem that the solvers can't do
anything about it. I assume this is due to the system matrix being singular. So
I transformed it to second order equations, which works fine in 1d (although
the system matrix is still singular. Weird...).
Now I went on to the 2d case and am now struggeling with the periodic boundary
conditions.
I'm using a hyper_rectangle grid and am applying Dirichlet conditions to the
left and right, i.e. face 0 and 1, and want to apply periodic boundary
conditions to the top and bottom, i.e. face 2 and 3.
To apply the periodic boundary conditions, I'm using the recursive approach
given in the step-45 tutorial. See the attachment for implementation.
In 1d the PDE can be solved analytically and in this specific configuration I
just expect the 1d solution periodically continued along the y-direction for
the E_x component and the E_y to vanish.
However there seems to be a problem with the periodic boundary conditions. When
multiplying the system matrix with the exact solution and subtracting the rhs
everything is fine except of the dofs on the constrained boundary where I
applied periodic boundary conditions, i.e. those on face 2.
Solving this matrix (using SolverGMRES) leads in deal.II to a solution for
which the periodic boundaries are fulfilled, but it's far away from the
expected solution.
I also exported the matrix and rhs to matlab and solved it there using gmres.
Here I get the expected solution apart from the dofs an face 2, which all have
the value 0. It looks like someone applied dirichlet boundaries there too.
I hope I could explain my problem understandable and now asking you for some
help on that to solve this problem.
If you need more information just ask, I will provide them.
Thanks.
Timo
void HSSolver2d::makeGridAndDofs()
{
std::cout << "Making grid and dofs..." << std::endl;
const dealii::Point<2> p1(-25e-4,-2.5e-4);
const dealii::Point<2> p2(25e-4,2.5e-4);
dealii::GridGenerator::hyper_rectangle(tria, p1, p2);
//before refining we set indicators for dirichlet and periodic boundary
conditions
dofHandler.begin_active()->face(0)->set_boundary_indicator(1);
dofHandler.begin_active()->face(1)->set_boundary_indicator(1);
//refine grid
tria.refine_global(2);
dofHandler.distribute_dofs(fe);
std::cout << "Number of active cells: "
<< tria.n_active_cells ()
<< std::endl
<< "Dofs: " << dofHandler.n_dofs ()
<< std::endl;
constraints.clear();
//periodic boundaries
applyPeriodicBoundaryConditions();
std::vector<bool> comps(2, false);
comps[0] = true;
//in 1d just set dirichlet
dealii::VectorTools::interpolate_boundary_values(dofHandler, 1,
dealii::ZeroFunction<2>(2),
constraints,
comps);
//try setting dirichlet without VectorTools
//applyDirichletBoundaryConditions(); works but same 'bug'
constraints.close();
//now create SparsityPattern
dealii::CompressedSparsityPattern cSparsity(dofHandler.n_dofs(),
dofHandler.n_dofs());
dealii::DoFTools::make_sparsity_pattern(dofHandler, cSparsity);
constraints.condense(cSparsity);
sPattern.copy_from(cSparsity);
systemMatrix.reinit(sPattern);
solution.reinit(dofHandler.n_dofs());
rhs.reinit(dofHandler.n_dofs());
//interpolate exact solution on dofs
dealii::VectorTools::interpolate(dofHandler, SolutionFunction<2>(),
solution);
}
void HSSolver2d::applyPeriodicBoundaryConditions()
{
applyPeriodicBoundaryConditions(dofHandler.begin(0)->face(2),
dofHandler.begin(0)->face(3));
}
void HSSolver2d::applyPeriodicBoundaryConditions(const
dealii::DoFHandler<2>::face_iterator &f1,
const dealii::DoFHandler<2>::face_iterator &f2)
{
if(f1->has_children())
{
for(unsigned int c = 0; c < f1->n_children(); ++c)
{
applyPeriodicBoundaryConditions(f1->child(c),
f2->child(c));
}
}
else
{
const unsigned int nDofsPerFace =
dofHandler.get_fe().dofs_per_face;
std::vector<unsigned int> localDofs1(nDofsPerFace);
f1->get_dof_indices(localDofs1);
std::vector<unsigned int> localDofs2(nDofsPerFace);
f2->get_dof_indices(localDofs2);
for(unsigned int i = 0; i < nDofsPerFace; ++i)
{
if(!constraints.is_constrained(localDofs1[i]) &&
!constraints.is_constrained(localDofs2[i]))
{
constraints.add_line(localDofs1[i]);
constraints.add_entry(localDofs1[i],
localDofs2[i], 1.0);
}
}
}
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
