Hi Konrad, I may add as a comment (without going into details): Your problem is - as you correctly describe above - a mixed Poisson problem. You essentially have two options to solve it and the way you choose influences the way you need to treat the boundary conditions: in primal form and in mixed form. In your formulation you have natural boundary conditions (Dirichlet) on Phi which means that you need to prescribe them in the variational form. This will lead to a boundary integral on the relevant part of the boundary and result in an additional term on the right-hand side. The boundary condition for E is a flux condition and in your mixed form this (Neumann) condition is essential. This means that you need to enforce them on the Raviart-Thomas space as a hard constraints. There are essentially two options for this:
1. If you have given a vector field that describes E on the relevant part of the boundary you can use the function *VectorTools::project_boundary_values_div_conforming. *This will project the given vector field onto its normal component (since this is the one you must describe with RT elements in H(div)) and use the projected value as boundary condition for the normal flux of E 2. If you only are given the scalar value of the normal flux of E then you may have to set the relevant dofs in the RT-space by hand. This you can do through constraints that way: AffineConstraints<double> constraints; constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, constraints); const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int dofs_per_face = fe.dofs_per_face; std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell); std::vector<types::global_dof_index> local_dof_face_indices (dofs_per_face); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { cell->get_dof_indices(local_dof_indices); for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n) { cell->face(face_n)->get_dof_indices(local_dof_face_indices); if (cell->at_boundary(face_n)) { if (cell->face(face_n)->boundary_id()==id_bottom_boundary) { for (unsigned int i = 0; i<dofs_per_face; ++i) { // see documentation for these functions constraints.add_line (local_dof_face_indices.at(i)); constraints.set_inhomogeneity ( local_dof_face_indices.at(i), value_of_dof); } } else // add some constraints on other boundaries. Here: zero-flux { for (unsigned int i = 0; i<dofs_per_face; ++i) { constraints.add_line (local_dof_face_indices.at(i)); } } } } ++cell_n; } constraints.close(); Once you use the DoFTools::make_sparsity_pattern function you need to make sure to keep the inhomogenities: DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, /*keep_constrained_dofs = */ true); The same in the assembly routine: cell->get_dof_indices(local_dof_indices); constraints.distribute_local_to_global(local_matrix, local_rhs_local, local_dof_indices, system_matrix, system_rhs, /* use inhomogeneities for rhs */ true); Keep in mind that the dofs for RT elements are edge moments in 2D and face moments in 3D (and not nodal values). Best, Konrad -- 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/48e77a20-69a8-4108-a5bc-d61adf3fba95%40googlegroups.com.