Dear deal.ii,
Please suggest me in implementing 'w.n = u.n' as a boundary condition where
'u' is Stokes solution with 'FESystem<dim> fe_stokes(FE_Q<dim>(degree+1),
dim, FE_Q<dim>(degree), 1)' and 'w' for mesh velocity with 'FESystem<dim>
fe_move(FE_Q<dim>(degree), dim)'. As I have taken 'degree=1', 'fe_move' has
unit support only at geometric vertices. I want to know if there is a
direct way to get outward unit normal at a boundary vertex which is average
of normals of faces that share this vertex(The faces are flat as I have not
assigned any MapingQ<dim>). If this is not possible/straight_forward, the
normal at a quadrature point on the face from which I visit this vertex the
first time is good enough.
I attempted implementing in the following way
MappingQGeneric<dim> mapping(fe.degree);
Quadrature<dim - 1> face_quadrature(femove.
get_unit_face_support_points());
FEFaceValues<dim> fe_face_values_stokes(mapping, fe,
face_quadrature, update_values | update_normal_vectors);
const unsigned int n_face_dofs = face_quadrature.size();
std::vector<Vector<double>> stokes_values(n_face_dofs, Vector
<double>(dim + 1));
std::vector<types::global_dof_index> local_dof_indices(femove.
dofs_per_face);
Tensor<1, dim> normal_vector1;
Tensor<1, dim> col_indices;
Tensor<1, dim> stokes_value_atvertex;
auto cell_move = move_dof_handler.begin_active();
auto endcell_move = move_dof_handler.end();
auto cell_stokes = dof_handler.begin_active();
std::vector<bool> vertex_touched(triangulation.n_vertices(),
false);
for (; cell_move != endcell_move; ++cell_stokes, ++cell_move)
{
if (!cell_move->is_artificial())
{
for (unsigned int f = 0; f < GeometryInfo<dim>::
faces_per_cell; ++f)
{
if (cell_move->face(f)->boundary_id() == 1006)
{
fe_face_values_stokes.reinit(cell_stokes, f);
normal_vector1 = fe_face_values_stokes.
normal_vector(1);
if(fabs(normal_vector1[0]) < 1e-10)
normal_vector1[0] = 1e-10;
fe_face_values_stokes.get_function_values(
lr_solution, stokes_values);
cell_move->face(f)->get_dof_indices(
local_dof_indices);
for (unsigned int v=0; v < GeometryInfo<dim>::
vertices_per_face; ++v)
{
if(vertex_touched[cell_move->face(f)->
vertex_index(v)] == false)
{
vertex_touched[cell_move->face(f)->
vertex_index(v)] = true;
for(int j=0; j < dim; ++j)
{
col_indices[j] = cell_move->face(f
)->vertex_dof_index(v, j);
if(j==0)
moveconstraints.add_line(
col_indices[0]);
const unsigned int component =
femove.system_to_component_index(v*dim+j).first;
stokes_value_atvertex[j] =
stokes_values[v*dim+j](component);
if(j!=0)
moveconstraints.add_entry(
col_indices[0], col_indices[j], normal_vector1[j]/normal_vector1[0]);
}
moveconstraints.set_inhomogeneity(
col_indices[0], -stokes_value_atvertex*normal_vector1/normal_vector1[0]);
}
}
}
}
}
}
but I get the following error.
void dealii::AffineConstraints<number>::add_entry(dealii::AffineConstraints
<number>::size_type, dealii::AffineConstraints<number>::size_type, number) [
with number = double; dealii::AffineConstraints<number>::size_type =
unsigned int]
The violated condition was:
std::abs(p.second - value) < 1.e-14
Additional information:
The entry for the indices 1143 and 1144 already exists, but the values
0.387236 (old) and 0 (new) differ by -0.387236.
Please correct any gaps in my understanding and suggest if there is a
general way which works with any order finite element.
Thanks,
Bhanu.
--
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/d9daac81-4d75-4724-883c-0ccd45391955%40googlegroups.com.
