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
