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
