Dear All,
I'm adapting tutorial step-27 to use discontinuous finite elements.
When assembling the local matrices to the system matrix step-27 uses
constraints.distribute_local_to_global (ui_vi_matrix, cell_rhs,
local_dof_indices,
system_matrix, system_rhs);
where ui_vi_matrix comes from assembling terms on the current cell
(the ui_vi notation being consistent with step-12).
Assuming that I want to continue to use
constraints.distribute_local_to_global, how do I treat the assembly
terms ui_ve_matrix, ue_vi_matrix i.e. those that contain integrals
involving cell and neighbor? The two cells may not have the same
number of degrees of freedom, so the matrices may not be square, so I
tried to use
void ConstraintMatrix::distribute_local_to_global ( const
FullMatrix< double > & local_matrix,
const std::vector< unsigned int > & row_indices,
const std::vector< unsigned int > & col_indices,
MatrixType & global_matrix
) const
but the local_dof_indices and neighbor_dof_indices are not the correct
variables to pass - how do I get row_indices and col_indices from the
dof_indices?
I have copied the assemble_system routine below for reference.
Thanks,
John
template <int dim>
void cdGMethod<dim>::assemble_system ()
{
// Define update flags for cells and faces
const UpdateFlags cell_update_flags = update_values
| update_gradients
| update_quadrature_points
| update_JxW_values;
const UpdateFlags face_update_flags = update_values
| update_quadrature_points
| update_JxW_values
| update_normal_vectors
| update_gradients;
// Neighbour normals and JxW values are same as on current cell, so
we don't
// update them
const UpdateFlags neighbor_face_update_flags = update_values
| update_gradients;
// create FEValues etc using default MappingQ1. Subface values not
needed
// currently.
hp::FEValues<dim> hp_fe_v (fe_collection,
quadrature_collection,
cell_update_flags);
hp::FEFaceValues<dim> hp_fe_v_face ( fe_collection,
face_quadrature_collection,
face_update_flags);
hp::FESubfaceValues<dim> hp_fe_v_subface ( fe_collection,
face_quadrature_collection,
face_update_flags);
hp::FEFaceValues<dim> hp_fe_v_face_neighbor(fe_collection,
face_quadrature_collection,
neighbor_face_update_flags);
// We cannot now define the size of the matrices a priori as we do
not know
// which finite element we use and this can change from cell to
cell. We
// must therefore resize the matrix at each step depending on the fe
used.
// i is internal, e external for shape functions u and test fns v.
FullMatrix<double> ui_vi_matrix;
FullMatrix<double> ue_vi_matrix;
FullMatrix<double> ui_ve_matrix;
FullMatrix<double> ue_ve_matrix;
Vector<double> cell_rhs;
std::vector<unsigned int> local_dof_indices;
std::vector<unsigned int> neighbor_dof_indices;
typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
// Loop over all cells
for (; cell!=endc; ++cell)
{
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
// Resize this cell matrix and rhs and reinit
ui_vi_matrix.reinit (dofs_per_cell, dofs_per_cell);
ui_vi_matrix = 0;
cell_rhs.reinit (dofs_per_cell);
cell_rhs = 0;
hp_fe_v.reinit (cell);
// Pick out the local fe_values
const FEValues<dim> &fe_v = hp_fe_v.get_present_fe_values ();
// Assemble the cell term
assemble_cell_term (fe_v, ui_vi_matrix, cell_rhs);
// Get the dof_indices for the current cell - we will use them later
local_dof_indices.resize (dofs_per_cell);
cell->get_dof_indices (local_dof_indices);
// Loop over all faces of current cell
for (unsigned int face_no=0;
face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
typename hp::DoFHandler<dim>::face_iterator
face=cell->face(face_no);
// Treat the boundary faces first - we only get one chance at
these
// so we must sum them when we find them
if (face->at_boundary())
{
hp_fe_v_face.reinit (cell, face_no);
// Pick out the local values
const FEFaceValues<dim> &fe_v_face
= hp_fe_v_face.get_present_fe_values ();
assemble_boundary_term (fe_v_face, ui_vi_matrix, cell_rhs);
}
else // internal face
{
Assert (cell->neighbor(face_no).state() == IteratorState::valid,
ExcInternalError());
typename hp::DoFHandler<dim>::cell_iterator
neighbor=cell->neighbor(face_no);
// Check we haven't done any refinement as this code hasn't been
// written yet
Assert (neighbor->level() == cell->level(), ExcNotImplemented());
// To make sure we only do each face once, only sum if
neighbors index
// is higher
if (neighbor->index() > cell->index())
{
// Determine the current cell in the neighbors world
const unsigned int
neighbor2=cell->neighbor_of_neighbor(face_no);
// Resize each matrix
const unsigned int neighbor_dofs_per_cell
=neighbor->get_fe().dofs_per_cell;
ui_ve_matrix.reinit (dofs_per_cell, neighbor_dofs_per_cell);
ui_ve_matrix = 0;
ue_vi_matrix.reinit (neighbor_dofs_per_cell, dofs_per_cell);
ue_vi_matrix = 0;
ue_ve_matrix.reinit (neighbor_dofs_per_cell,
neighbor_dofs_per_cell);
ue_ve_matrix = 0;
hp_fe_v_face.reinit (cell, face_no);
hp_fe_v_face_neighbor.reinit (neighbor, neighbor2);
// Pick out the local values
const FEFaceValues<dim> &fe_v_face
= hp_fe_v_face.get_present_fe_values ();
const FEFaceValues<dim> &fe_v_face_neighbor
= hp_fe_v_face_neighbor.get_present_fe_values ();
assemble_face_term (fe_v_face,
fe_v_face_neighbor,
ui_vi_matrix,
ue_vi_matrix,
ui_ve_matrix,
ue_ve_matrix);
// Write to the system matrix for things involving the neighbor
neighbor_dof_indices.resize (neighbor_dofs_per_cell);
neighbor->get_dof_indices (neighbor_dof_indices);
// Want to distribute the local matrix to the global, but
these calls
// are not correct
constraints.distribute_local_to_global (ui_ve_matrix,
local_dof_indices,
neighbor_dof_indices,
system_matrix);
constraints.distribute_local_to_global (ue_vi_matrix,
neighbor_dof_indices,
local_dof_indices,
system_matrix);
// This should be fine, though
constraints.distribute_local_to_global (ue_ve_matrix,
neighbor_dof_indices,
system_matrix);
}
}
} //end of loop over faces
// Now write ui_vi_matrix to the system matrix, including rhs in
the same
// way as for the continuous version
constraints.distribute_local_to_global (ui_vi_matrix, cell_rhs,
local_dof_indices,
system_matrix, system_rhs);
} //end of loop over cells
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
