Your suggestions worked exactly as advertised. I've copy-pasta's my
updated code below for completeness.
> 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.
>
>
template <int spacedim>
void VectorHelfrichFlow<spacedim>::data_dot_normal (const Vector<double>
&data)
{
> > {
> > /*{{{*/
> > 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
>
> std::vector<Tensor<1,spacedim> > local_data_values(n_q_points,
Tensor<1,spacedim>());
const FEValuesExtractors::Vector W (0);
unsigned int k=0;
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[W].get_function_values(data,local_data_values);
double data_dot_normal_at_q_points = 0;
for (unsigned int q=0; q<n_q_points; ++q)
}
data_dot_normal_at_q_points +=
local_data_values[q]*fe_values.normal_vector(q);
}
double data_dot_normal_on_cell = data_dot_normal_at_q_points/n_q_points;
> > 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.
>
>
cellwise_H(k) = data_dot_normal_on_cell;
k+=1; // k enumerates the loop through cells
Thank you for your help.
--
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.
