Thomas,
your code looks essentially correct, though it can be improved:
Now, for my attempt at this:
|
template(int spacedim>
void data_dot_normal (Vector<double> data)
This copies the entire 'data' vector. You'll probably be better off making it
a 'const' reference.
{
/*{{{*/
MappingQEulerian<dim,Vector<double>,spacedim> mapping(2, dof_handler,
euler_vector);
H = 0;
const QGauss<dim> quadrature_formula (2*fe.degree);
FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula,
update_normal_vectors);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<Vector<double> > local_data_values(n_q_points,
Vector<double>(spacedim));
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
cell!=endc; ++cell)
{
fe_values.reinit(cell);
fe_values.get_function_values(data,local_data_values);
You can improve this a bit if you create an FEValuesExtractor::Vector, and use
this here to get 'local_data_values' not as a std::vector<Vector<double>>,but
as a std::vector<Tensor<1,spacedim>. If you have this, then you can also...
double data_dot_normal_at_q_points = 0;
for (unsigned int q=0; q<n_q_points; ++q)
{
Point<spacedim> space_point = fe_values.quadrature_point(q);
Tensor<1,spacedim> unit_normal = fe_values.normal_vector(q);
double summ = 0;
for (unsigned int d=0; d<spacedim; ++d)
{
summ += local_data_values[q][d]*unit_normal[d];
}
...replace this loop here by simply saying
local_data_values[q] * unit_normal
data_dot_normal_at_q_points += summ;
}
double data_dot_normal_on_cell = data_dot_normal_at_q_points/n_q_points;
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
H(local_dof_indices[i]) = data_dot_normal_on_cell;
}
}
}
|
For each cell, I am dotting the vector-valued data with the normal at each
quadrature point, and since the result of that is only defined at dofs, I
average the n_q_points-many dot products and fill the result with that average.
There are two things wrong with this:
1. Interestingly, this appears essentially correct when plotted, but for
coarse meshes the scalar field is skewed over my surface ever so slightly,
and for finer meshes the scalar field begins to line up with where I
expect it to on my geometry. I'm not sure what that's about.
This is because you compute a single scalar number for each cell, but then you
write it into the vector locations for *all* degrees of freedom associated
with this cell. If you visit the next cell, you're going to overwrite the
values you previously put into some of the degrees of freedom. Which of the
cells adjacent to a vertex wins depends on the order of cells. In your case,
they seem to be ordered in a particular direction, leading to what you
describe as a "skew".
2. The (scalar) result, H, is defined at each dof (not surprising), but I
don't need it at each dof, in fact H has three times as many dofs as it
needs. How can I specify a scalar field over my surface along was well as
the vector-valued computed solution?
You can declare a
Vector<double> cellwise_H (triangulation.n_active_cells());
and then fill that in the order in which you traverse the cells, one element
per cell. You can use that for further computations. You can even pass this
vector to DataOut.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.colostate.edu/~bangerth/
--
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.
