Here's one way around this:
At the end of the day, the problem I had is that I did not understand how a
Manifold object interacts with a Triangulation and a Mapping object. I
probably still don't, so take this with a grain of salt:
I was creating a Manifold object and then attaching it to a Triangulation,
like so
static Ellipsoid<dim,spacedim> ellipsoid(a,b,c,center);
static SphericalManifold<dim,spacedim> sphere();
GridGenerator::hyper_sphere(triangulation,center, 1);
triangulation.set_all_manifold_ids(0);
triangulation.set_manifold (0, sphere);
After setting up the DoFHandler and FESystem, a call to
FEValues<dim,spacedim> fe_values(mapping, fe, quadrature_formula,...);
maps all points from the triangulation -- *and quadrature points* --
according to the push_forward function in SphericalManifold. The result is
that we can easily compute shape gradients on this surface.
Alternatively, when the Mapping is not MappingQ, but MappingQEulerian, the
mapping waits for specific instructions on how to move the quadrature
points to the desired shape - that sits in the euler_vector. I was passing
a zero euler_vector because I figured that my manifold and triangulation
would have figured out what to do with all relevant points.
The fix here is to get rid of the SphericalManifold altogether, and just
start with a triangulation. Then create an euler_vector by interpolating a
function which can map arbitrary points to the desired geometry (this is
essentially the push_forward function of SphericalManifold if you want a
sphere). This euler_vector will not be zero!
The whole thing looks like this:
GridGenerator::hyper_sphere(triangulation,center,1.0);
euler_vector.reinit (dof_handler.n_dofs());
VectorTools::interpolate(dof_handler,
MapToSphere<spacedim>(r),
euler_vector);
MappingQEulerian<dim,Vector<double>,spacedim> mapping(2, dof_handler,
euler_vector);
FEValues<dim,spacedim> fe_values (mapping, fe, quadrature_formula, ...)
Now one can obtain a sphere, compute shape derivatives on it, and then
update the euler_vector in order to deform the sphere into something else.
I suggest to myself to reread the Interplay of UpdateFlags, Mapping, and
FiniteElement in FEValues:
https://dealii.org/8.4.1/doxygen/deal.II/group__UpdateFlags.html, maybe
this will all become clear one day ;)
On Wednesday, October 19, 2016 at 7:19:47 PM UTC-4, thomas stephens wrote:
>
> I am having trouble understanding how to use MappingQEulerian - I have
> constructed a codimension-1 sphere by attaching a SphericalManifold to a
> hyper_sphere triangulation.
>
> So far so good.
>
> Now, I have some surface derivatives I have been computing within the
> function below called assemble_system()
>
> assemble_system ()
> {
> /*{{{*/
>
> /*
> * [ M L-2*hL+0.5*d ][ V_n+1 ] [ 0 ]
> * | || | = | |
> * [ -zn*L M ][ H_n+1 ] [ rhs ]
> *
> * system_matrix*VH = system_rhs
> *
> * system_matrix has size 2*n_dofs x 2*n_dofs,
> * each block has size n_dofs x n_dofs, and
> * rhs has size n_dofs
> *
> */
>
>
> *MappingQ<dim,spacedim> mapping(2);*
>
> system_matrix = 0;
> VH = 0;
> system_rhs = 0;
>
> 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);
>
> and life has been pretty good. The subsequent derivatives look like this:
>
> const *FEValuesExtractors::Vector W (0)*;
> std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
>
> for (typename DoFHandler<dim,spacedim>::active_cell_iterator
> cell = dof_handler.begin_active(),
> endc = dof_handler.end();
> cell!=endc; ++cell)
> {
> local_M = 0;
> local_L = 0;
> local_hL = 0;
> local_d = 0;
> local_rhs = 0;
>
> fe_values.reinit (cell);
>
> for (unsigned int i=0; i<dofs_per_cell; ++i)
> for (unsigned int j=0; j<dofs_per_cell; ++j)
> for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
> {
>
> local_M(i,j) +=
>
> *fe_values[W].value(i,q_point)*
> fe_values[W].value(j,q_point)*
> fe_values.JxW(q_point)*;
>
> local_L(i,j) +=
>
> *scalar_product(fe_values[W].gradient(i,q_point),
>
> fe_values[W].gradient(j,q_point) )*
> fe_values.JxW(q_point)*;
>
> local_hL(i,j) +=
>
>
> *scalar_product(2.0*identity_on_manifold.shape_grad(fe_values.normal_vector(q_point))*
>
> fe_values[W].gradient(i,q_point),
> fe_values[W].gradient(j,q_point) )*
> fe_values.JxW(q_point)*;
>
> But now I want to move my mesh, and my strategy has been to compute a
> deformation vector (of size dof_handler.n_dofs()), and create a new
> MappingQEulerian at each time step using this updated deformation vector.
>
> My problem is that I cannot really get off the ground with this strategy
> because (it seems to me that) once I replace MappingQ with
> MappingQEulerian, these derivatives are no longer correct. So, for the
> modified assemble_system() code I have:
>
> assemble_system ()
> {
> /*{{{*/
>
> /*
> * [ M L-2*hL+0.5*d ][ V_n+1 ] [ 0 ]
> * | || | = | |
> * [ -zn*L M ][ H_n+1 ] [ rhs ]
> *
> * system_matrix*VH = system_rhs
> *
> * system_matrix has size 2*n_dofs x 2*n_dofs,
> * each block has size n_dofs x n_dofs, and
> * rhs has size n_dofs
> *
> */
>
> /* build a MappingQEulerian that approximates
> * grid_transform : sphere --> ellipsoid. This will be updated by the
> * deformation induced by the velocity field computed at each time step.
> */
>
> /* create an interpolation of the transformation that takes the sphere
> to the
> * ellipsoid, this will be updated at each mesh movement */
>
> deformation = 0;
> * MappingQEulerian<dim,Vector<double>,spacedim> mapping(2, dof_handler,
> deformation);*
>
>
> system_matrix = 0;
> VH = 0;
> system_rhs = 0;
>
> 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);
>
>
> The difference in output looks like this:
>
>
> From *MappingQ*:
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> From *MappingQEulerian*:
>
>
>
>
>
>
> Something is very wrong here - I specifically set deformation=0; before
> building MappingQEulerian in order to test the null deformation (i.e. the
> identity), but I still get some silly numbers for H2 (which is computed
> directly from those derivatives.
>
> What am I missing? Let me know if there is more information I can provide
> here.
>
>
> Thank you,
> 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.
