Well, here's a pretty ugly attempt -
One thing I want to point out is that my solution vector comes from a
finite element built like this: (dim=2, spacedim = 3)
FESystem<dim,spacedim> fe;
fe(FE_Q<dim,spacedim>(fe_degree), spacedim)
and the only DoFHandler I have lying around right now is
dof_handler.distribute_dofs(fe)
so my dofs correspond to a 3-component element on a surface.
Furthermore, my surface is defined through an euler_vector passed to
MappingQEulerian, so computations using normal vectors should use the
deformed normal vectors - obtained through something like this:
MappingQEulerian<dim,Vector<double>,spacedim> mapping(2, dof_handler,
euler_vector);
const QGauss<dim> quadrature_formula (2*fe.degree);
FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula,
update_values |
update_normal_vectors |
update_gradients |
update_quadrature_points |
update_JxW_values);
...
fe_values.normal_vector(q_point); //q_point being the index into some
current cell's list of quadrature points
*My goal* is to use my vector-valued solution to compute a scalar quantity
on my surface which depends on normal vectors, and ultimately make some
computational decisions with the result. (I would also like to plot the
scalar result of my computation.)
Now, for my attempt at this:
template(int spacedim>
void data_dot_normal (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);
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];
}
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.
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?
To Wolfgang, perhaps I should go back to DataPostProcessor, and perhaps I
should try to activate those normals, let me know what you think.
Tom
On Wednesday, October 26, 2016 at 5:30:51 PM UTC-4, thomas stephens wrote:
>
> I have a vector-valued solution, call it H, on a codimension-one geometry
> (described by a MappingQEulerian object) and I am very interested in
> obtaining H*nu, the dot product of H with the normal vector nu.
>
> I have already obtained H*H using a class derived from
> DataPostProcessorScalar, and that makes sense to me - the function
> signature of compute_derived_quantities_vector() looks like this:
>
> template <int spacedim>
> void
> VectorValuedSolutionSquared<spacedim>::compute_derived_quantities_vector
> (const std::vector<Vector<double> > &uh,
> const
> std::vector<std::vector<Tensor<1, spacedim> > > & /*duh*/,
> const
> std::vector<std::vector<Tensor<2, spacedim> > > & /*dduh*/,
> const
> std::vector<Point<spacedim> > & /*normals*/,
> const
> std::vector<Point<spacedim> > & /*evaluation_points*/,
>
> std::vector<Vector<double> >
> &computed_quantities) const
>
>
> since H*H doesn't require anything other than uh, but now I would like to
> take advantage of that vector of normals (which, for dealii 8.44.1,
> should be of type const std::vector<Tensor<1,spacedim> >, see here - but I
> digress)
>
> At any rate, I don't know how to access this normals vector, and if I
> were to compute it ahead of time using fe_values.get_all_normals() then I
> am not sure how to pass it to compute_derived_quantities(), since my code
> never actually calls that function.
>
> At the end of the day I think this comes down to me fighting with a
> solution H of type Vector<double> of size dof_handler.n_dofs(), where that
> dof_handler was built with vector-valued finite elements, and now I want to
> look at essentially a scalar function on my geometry.
>
> Any help would be great,
> Tom
>
>
--
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.
