Hello all,
I am using deal.ii to solve two vector valued equations by FEM:- one for
the velocity field "*v*" and the other for a vector field "*f*".
I created different classes for solving these two equations separately.
I am using same parallel::distributed::triangulation for both the equations
but they have different DoFHandler object.
After solving for the velocity field "v", I need to multiply the stiffness
of Velocity equation to its velocity field "v" to get the global reactions
= Rv.
For the second equation, the "global Rhs vector" comes from two parts:
1) R1 : part one comes from the weak form of the equation for "f" .
2) part two is Rv which is calculated above from the solution of velocity
equation.
So, if the solution variable is f, the final discretized version of
equation 2 becomes
RHS = R = R1 + Rv
[K] {f} = R1 + Rv
where K is the system stiffness matrix for f equation.
Following examples, I assemble the stiffness and "R1" using the following
commands
constraints.clear();
constraints.reinit (locally_relevant_dofs);
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
// constraints do not include anything else
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
R1);
To get R, R1 and Rv should be added in such a way that at any vertex of the
triangulation,
if the dofs from the velocity equation are dv1, dv2, dv3 and dofs from the
"f" equation are df1, df2, df3
then
R(df1) = R1(df1) + Rv(dv1)
R(df2) = R1(df2) + Rv(dv2)
R(df3) = R1(df3) + Rv(dv3)
Here, I am assuming that since the dof handler object for two equations are
different, the dof number at a vertex could be different.
I am unable to come out with a way to assemble the Rhs vector "R" which is
sum of two global vectors, "R1" and "Rv".
I have been unable to figure out how to do this for quite a while now, so
any help will be appreciated.
Here is what I have come up until now,
but there are issues which I am unable to figure out how to solve, like,
Dofs that are at the junction of two processors are accounted for twice
which is wrong.
typename DoFHandler<dim>::active_cell_iterator
fcell = fdof_handler.begin_active(),
fendc = fdof_handler.end(),
v_cell = veqn->dof_handler.begin_active();
std::map<unsigned int, bool> vertex_touched;
for (; fcell != fendc; ++fcell)
if (fcell->is_locally_owned())
for (unsigned int v = 0; v < GeometryInfo<dim>::
vertices_per_cell; ++v)
vertex_touched[fcell->vertex_index(v)] = false;
fcell = dof_handler.begin_active();
for(; fcell!= fendc; ++v_cell,++cell)
if(v_cell->is_locally_owned())
if(v_cell->at_boundary())
for(unsigned int face = 0; face < GeometryInfo<dim>::
faces_per_cell; ++face)
{
if(v_cell->face(face)->at_boundary())
for(unsigned int v = 0; v < GeometryInfo<dim>::
vertices_per_face; ++v)
{
if(vertex_touched[fcell->vertex_index(v)] == false)
{
vertex_touched[fcell->vertex_index(v)] == true;
for(unsigned int d = 0; d < dim; ++d)
{ R(fcell->face(face)->vertex_dof_index(v,d))
= R1(fcell->face(face)->vertex_dof_index(v,d)) +
veqn->Rv(v_cell->face(face)->
vertex_dof_index(v,d));
}
}
}
}
--
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].
For more options, visit https://groups.google.com/d/optout.
