Hello all,
I have a problem related to step-33 example, which I am modifying. This
example uses a shock capturing term, which I am now modifying according to
the thesis of Hartmann. The term I want is of the form
integral( eps * grad(u) * grad(v) )
with v being the test function, where
eps = C * h^p * | div(F(u)) |
with F being the flux matrix. I compute the flux divergence as follows. The
shocks were not coming out correctly, so I tested to see if the divergence
was non-zero, and it turns out that it is close to machine zero. So there is
some mistake in this code but I cannot find it. Hope somebody can help me.
F is a vector of size EulerEquations<dim>::n_components
Thanks
praveen
Table<2,Sacado::Fad::DFad<double> >
divergence(n_q_points, EulerEquations<dim>::n_components);
// Compute divergence of flux
for (unsigned int point=0; point<n_q_points; ++point)
{
for (unsigned int c=0; c<EulerEquations<dim>::n_components; ++c)
divergence[point][c] = 0.0;
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int
component_i = fe_v.get_fe().system_to_component_index(i).first;
for(unsigned int d=0; d<dim; ++d)
divergence[point][component_i] += flux[point][component_i][d] *
fe_v.shape_grad_component(i,
point, component_i)[d];
}
for (unsigned int c=0; c<EulerEquations<dim>::n_components; ++c)
if(std::fabs(divergence[point][c]) > 1.0e-10)
std::cout << "div = " << (divergence[point][c]) << "\n";
}
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
