I found the mistake in my code. Please ignore this issue. Thanks praveen
On Sat, Apr 2, 2011 at 8:32 PM, Praveen C <[email protected]> wrote: > 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
