Hello
My fe space is
FE_DGQArbitraryNodes<dim> (QGauss<1>(degree+1))
For face quadrature I use
QGauss<dim-1>(degree+1)
At each face quadrature point, there are only (degree+1) dofs which are
supported there. I would like to find this information. I have written a
function as below which computes
Table<3,unsigned int> face_dof_index_u;
face_dof_index_u[f][q] gives all dofs supported on f'th face and q'th
quadrature point.
Is such information available in the library ?
I am currently using following function to find this information. Is there
a better way to find it ?
This information helps me to write more efficient face assembly. I also
need to find Jacobian for Newton method which can be made more efficient
with such information.
Thanks
praveen
//------------------------------------------------------------------------------
// Find which dofs have non-zero support on each quadrature on faces
//------------------------------------------------------------------------------
template <int dim>
void
ConservationLaw<dim>::find_face_dof_indices ()
{
pcout << "Finding face dof indices\n";
TableIndices<3> size_u (GeometryInfo<dim>::faces_per_cell,
fe_u.degree+1,
fe_u.degree+1);
face_dof_index_u.reinit (size_u);
// Tolerance to decide support of basis function
const double TOL = 1.0e-12;
QGauss<dim-1> quadrature_formula (fe_u.degree+1);
const unsigned int n_q_points = quadrature_formula.size();
FEFaceValues<dim> fe_face_values_u (mapping(), fe_u,
quadrature_formula,
update_values);
// use the first cell
typename DoFHandler<dim>::active_cell_iterator
cell_u = dof_handler_u.begin_active();
for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
fe_face_values_u.reinit (cell_u, f);
for(unsigned int q=0; q<n_q_points; ++q)
{
unsigned int nu = 0;
for(unsigned int i=0; i<fe_u.dofs_per_cell; ++i)
if(fabs(fe_face_values_u.shape_value(i, q)) > TOL)
{
face_dof_index_u[f][q][nu] = i;
++nu;
}
AssertThrow(nu == fe_u.degree+1,
ExcMessage("Wrong no. of face U dofs"));
}
}
}
--
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.
