Dear all,
The file attached contains functions that might be a solution to the
problem posted. They work well for me. I hope this would be helpful for
people who have not solved the problem.
Best regards,
Jie
On Monday, March 18, 2019 at 6:18:54 PM UTC+1 Wolfgang Bangerth wrote:
> On 3/16/19 7:12 AM, Konrad wrote:
> >
> > Tried it with BDM elements but it does not work. I strongly think about
> > starting to contribute to deal.ii.
>
> I'm glad to see that you've already found the issue, but I would like to
> add
> that I think this here would be a really good idea! :-) We are always
> happy to
> have new contributors!
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
>
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/734220fa-dd49-483c-87b7-df8fba720f88n%40googlegroups.com.
namespace internals
{
template <typename cell_iterator>
void
compute_face_projection_div_conforming_bdm
(const cell_iterator &cell,
const unsigned int face,
const FEFaceValues<2> &fe_values,
const unsigned int first_vector_component,
const Function<2> &boundary_function,
const DerivativeForm<1,2,2>
&jacobian_of_one_quadrature_point,
AffineConstraints<double> &constraints)
{
const FEValuesExtractors::Vector vec (first_vector_component);
const FiniteElement<2> &fe = cell->get_fe ();
const unsigned int
face_coordinate_direction[GeometryInfo<2>::faces_per_cell] = {0, 0, 1, 1};
std::vector<Vector<double> > values (fe_values.n_quadrature_points,
Vector<double> (2));
Vector<double> dof_values (fe.dofs_per_face);
{
const std::vector<Point<2> > &quadrature_points =
fe_values.get_quadrature_points ();
boundary_function.vector_value_list (quadrature_points, values);
}
double param_jacobian = 0.0;
param_jacobian = std::sqrt
(jacobian_of_one_quadrature_point[0][face_coordinate_direction[face]]
*
jacobian_of_one_quadrature_point[0][face_coordinate_direction[face]]
+
jacobian_of_one_quadrature_point[1][face_coordinate_direction[face]]
*
jacobian_of_one_quadrature_point[1][face_coordinate_direction[face]]);
for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
{
dof_values (i) =
values[i][face_coordinate_direction[face]]*param_jacobian;
}
std::vector<types::global_dof_index> face_dof_indices (fe.dofs_per_face);
cell->face (face)->get_dof_indices (face_dof_indices,
cell->active_fe_index ());
// Copy the computed values in the ConstraintMatrix only, if the degree
// of freedom is not already constrained.
for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
if (!(constraints.is_constrained (face_dof_indices[i])))
{
constraints.add_line (face_dof_indices[i]);
if (std::abs (dof_values (i)) > 1e-14)
{
constraints.set_inhomogeneity (face_dof_indices[i], dof_values
(i));
}
}
}
}
template <int dim>
void
project_boundary_values_div_conforming_bdm (const DoFHandler<dim>
&dof_handler,
const unsigned int
first_vector_component,
const Function<dim>
&boundary_function,
const types::boundary_id
boundary_component,
AffineConstraints<double>
&constraints,
const Mapping<dim> &mapping)
{
const unsigned int spacedim = dim;
const FiniteElement<dim> &fe = dof_handler.get_fe ();
QGauss<dim - 1> face_quadrature (fe.degree);
FEFaceValues<dim> fe_face_values (mapping, fe, face_quadrature,
update_JxW_values |
update_normal_vectors |
update_quadrature_points |
update_values);
hp::FECollection<dim> fe_collection (fe);
hp::MappingCollection<dim> mapping_collection (mapping);
hp::QCollection<dim> quadrature_collection;
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
quadrature_collection.push_back (QProjector<dim>::project_to_face
(face_quadrature,
face));
hp::FEValues<dim> fe_values (mapping_collection, fe_collection,
quadrature_collection,
update_quadrature_points | update_jacobians);
for (typename DoFHandler<dim>::active_cell_iterator cell =
dof_handler.begin_active (); cell != dof_handler.end (); ++cell)
{
if (cell->at_boundary ())
{
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face)
{
if (cell->face (face)->boundary_id () == boundary_component)
{
fe_values.reinit (cell, face + cell->active_fe_index () *
GeometryInfo<dim>::faces_per_cell);
const DerivativeForm<1,dim,spacedim>&
jacobian_of_one_quadrature_point = fe_values.get_present_fe_values
().get_jacobians ()[0];
fe_face_values.reinit (cell, face);
internals::compute_face_projection_div_conforming_bdm (cell,
face,
fe_face_values,
first_vector_component,
boundary_function,
jacobian_of_one_quadrature_point,
constraints);
}
}
}
}
}
