Hi,
I am a complete beginner and trying to solve the 1D advection equation. I am
using DG method with an upwind flux and was wondering about the following
questions when I was writing the code.
1. Does the package use legendre or lagrange polynomials or is there some
way to define which one to use?
2. I am trying to get the value at the end point of the neighbor element for
the flux. Right now I am trying to use Gauss Lobatto points and using the
last quadrature point to get the boundary value. Is there a simpler way to
access the values at the end points for each element?
3. I am confused with the "reinit" function. Does it matter if I don't give
this command. I wanted to make sure because every program seems to use this.
4. Below is the code I've written for assembling the system. When I try to
access get_function_values(), I get the folowing error:
Error message:
An error occurred in line <1889> of file <source/fe/fe_values.cc> in
function
void dealii::FEValuesBase<dim, spacedim>::get_function_values(const
InputVector&, std::vector<dealii::Vector<number3>,
std::allocator<dealii::Vector<number3> > >&) const [with InputVector =
dealii::BlockVector<double>, number = double, int dim = 1, int spacedim = 1]
The violated condition was:
fe_function.size() == present_cell->n_dofs_for_dof_handler()
The name and call sequence of the exception was:
ExcDimensionMismatch(fe_function.size(),
present_cell->n_dofs_for_dof_handler())
Additional Information:
Dimension 0 not equal to 32
Stacktrace:
-----------
#0 /home/Ganesh/Research/Software/deal.II/lib/libdeal_II_1d.g.so: void
dealii::FEValuesBase<1, 1>::get_function_values<dealii::BlockVector<double>,
double>(dealii::BlockVector<double> const&,
std::vector<dealii::Vector<double>, std::allocator<dealii::Vector<double> >
>&) const
#1 ./project_v2: Magma_Dynamics<1>::assemble_system()
#2 ./project_v2: Magma_Dynamics<1>::run()
#3 ./project_v2: main
--------------------------------------------------------
make: *** [run] Aborted
Note: 32 is the number of active cells for this code.
Code:
template <int dim>
void Magma_Dynamics<dim>::assemble_system()
{
QGaussLobatto<dim> quadrature_formula(4);
// QGauss<dim-1> face_quadrature_formula(2);
FEValues<dim> fe_values(fe,quadrature_formula,
update_values | update_gradients
| update_JxW_values);
FEValues<dim> neighbor_fe_values(fe,quadrature_formula,
update_values
|update_gradients
|
update_JxW_values);
// FEFaceValues<dim> fe_face_values(fe,face_quadrature_formula,
//
update_values|update_quadrature_points |
//
update_JxW_values | update_normal_vectors);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell,dofs_per_cell);
FullMatrix<double> stiff_matrix(dofs_per_cell,dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
Vector<double> face_values(dofs_per_cell);
std::vector<unsigned int> local_dof_indices(dofs_per_cell);
std::vector<Vector<double> >
old_solution_values(n_q_points,Vector<double>(1));
std::vector<Vector<double> >
old_solution_values_neighbor(n_q_points,Vector<double>(1));
std::vector<Vector<double> >
present_solution_values(n_q_points,Vector<double>(1));
std::vector<Vector<double> >
present_solution_values_face(n_q_points,Vector<double>(1));
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
// Normal defined only for 1 dimension, to be corrected if changing
dimensions
const Point<dim> normal = Point<dim>(1.0);
for(;cell!=endc;++cell)
{
fe_values.reinit(cell);
cell_matrix=0.;
cell_rhs=0.;
fe_values.get_function_values(old_solution,old_solution_values);
fe_values.get_function_values(solution,present_solution_values);
for(unsigned int i=0;i<dofs_per_cell;++i)
for(unsigned int j=0;j<dofs_per_cell;++j)
{
for(unsigned int q_point=0;q_point<n_q_points;++q_point)
cell_matrix(i,j) += (fe_values.shape_value(i,q_point)*
fe_values.shape_value(j,q_point)*
fe_values.JxW(q_point));
for(unsigned int q_point=0;q_point<n_q_points;++q_point)
stiff_matrix(i,j) = (normal*fe_values.shape_value(i,q_point)*
fe_values.shape_grad(j,q_point)*
fe_values.JxW(q_point));
}
if (cell->at_boundary(0) == false)
{
typename DoFHandler<dim>::cell_iterator
left_neighbor= cell->neighbor(0);
while (left_neighbor->has_children())
left_neighbor = left_neighbor->child(1);
neighbor_fe_values.reinit (left_neighbor);
neighbor_fe_values.get_function_values(old_solution,
old_solution_values_neighbor);
}
}
}
Thanks a lot.
Ganesh
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
