There is nothing wrong with what you are doing.
The problem is in the nature of your Manifold, since it contains 3 singular
points. The first is the center of the ellipsoid, while the second and third
are the north and south poles.
Try attaching a SphericalManifold to your deformed grid, and see if you like
the result. If you do, then hack the SphericalManifold. I don’t know what
version of deal.II you are using. Assuming you are not using the latest dev, if
you open up the definition of spherical manifold in manifold_lib.cc, you’ll see
that SphericalManifold::get_new_point does something special in the case
spacedim == 3, i.e., SphericalManifold is not *really* a ChartManifold in 3d.
In this case the ChartManifold mechanism cannot be used for the reason above.
A quick explanation why things are not working:
Consider the top cell, and let us assume for a second that the north pole
(z=z_max) is in the center of the cell, and that the four vertices of your
ellipsoid lie at the same z (z=h < z_max). When you transform using
ChartManifold, the z gets mapped to the angle phi, and every point of the cell
will have the same angle phi, and different theta. If you take the average of
the four points in the chart manifold coordinates, you’ll see that the average
will have an average theta (which is ok), but it will also have the same phi of
the four surrounding points… in other words, the average will not be in the
middle of the four points (that should be the north pole: a singularity of your
mapping). It will lye on the same latitude of the other four points.
ChartManifold is doing exactly what you asked it to do, but in the case of maps
with singularities, you should not use ChartManifold…
My suggestion is to derive EllipsoidManifold from SphericalManifold, and then
overload directly get_new_point, by calling internally get_new_point of
spherical manifold, and then projecting the result to the ellipsoid...
Best,
Luca.
> On 09 Aug 2016, at 24:13, thomas stephens <tdstephe...@gmail.com> wrote:
>
>
>
> I am trying to obtain the mesh for a codimension-1 ellipsoid and attach an
> ellipsoidal manifold to it in order to refine the mesh. My strategy is
> failing and I have a few questions.
>
> My strategy:
> some definitions:
> dim=2; spacedim=3; chartdim=spacedim-1; a=1; b=2;c=3;
>
> Set up codimension-1 Triangulation:
>
> Triangulation<dim,spacedim> tria;
>
>
> Use GridGenerator::hyper_sphere() to obtain a codimension 1 mesh:
>
>
> GridGenerator::hyper_sphere (tria, Point<spacedim>(0.0,0.0,0.0), 1.0);
>
> Use GridTools::transform() to map triangulation of sphere onto triangulation
> of ellipsoid (I have a question about this, see below). Here, grid_transform
> has signature Point<spacedim> grid_transform(const Point<spacedim> X)
>
> GridTools::transform(&grid_transform, tria);
>
> Next, subclass ChartManifold<dim,spacedim,chartdim> in order to obtain a
> push_forward() and a pull_back() function for my parametric ellipsoidal
> manifold description. Then attach that manifold to the triangulation:
>
> Ellipsoid<dim,spacedim,chartdim> ellipsoid(a,b,c);
> tria.set_manifold(0,ellipsoid);
>
> Now refine the triangulation in order to see a refined ellipsoidal mesh:
>
> tria.refine_global(1);
>
> The Result: Complete garbage! I also don't really have an idea why this is
> not working - I've checked the dimensions on my push_forward() and
> pull_back() functions, I think that I'm setting the periodicity correctly in
> my ChartManifold superclass, and I've also checked that my simple math is
> okay. Let me know what other information would be useful. Below is the
> refined mesh for a=b=c=1. This should just be a sphere.
>
>
>
>
> The code:
> template <int dim,int spacedim,int chartdim=spacedim-1>
> class Ellipsoid: public ChartManifold<dim,spacedim,chartdim>
> {
> public:
>
> Ellipsoid(double,double,double);
>
> Point<chartdim> pull_back(const Point<spacedim> &space_point) const;
>
> Point<spacedim> push_forward(const Point<chartdim> &chart_point) const;
>
> private:
> double a,b,c;
> double max_axis;
> const Point<spacedim> center;
> SphericalManifold<dim,spacedim> reference_sphere;
>
> };
>
>
> template <int dim, int spacedim, int chartdim>
> Ellipsoid<dim,spacedim,chartdim>::Ellipsoid(double a, double b, double c):
> ChartManifold<dim,spacedim,chartdim>(Point<chartdim>(2*numbers::PI,
> 2*numbers::PI)), a(a), b(b),c(c), center(0,0,0), reference_sphere(center)
> {
> max_axis = std::max(std::max(a,b),c);
> }
> template <int dim,int spacedim, int chartdim>
> Point<chartdim> Ellipsoid<dim,spacedim,chartdim>::pull_back(const
> Point<spacedim> &space_point) const
> {
> double x,y,z, u,v,w;
>
> // get point on ellipsoid
> x = space_point[0];
> y = space_point[1];
> z = space_point[2];
>
> std::cout << "using a,b,c: " << std::endl;
> std::cout << a << " " << b << " " << c << std::endl;
> std::cout << "from pull_back: " << std::endl;
> std::cout << "space_point: " << std::endl;
> std::cout << x << " " << y << " " << z << std::endl;
>
> // map ellipsoid point onto sphere
> u = x/a;
> v = y/b;
> w = z/c;
>
> std::cout << "pulls back to : " << std::endl;
> std::cout << u << " " << v << " " << w << std::endl;
> std::cout << "on sphere." << std::endl;
>
> Point<spacedim> p(u,v,w);
>
> // use reference_sphere's pull_back function
> Point<spacedim> q = reference_sphere.pull_back(p);
> Point<chartdim> chart_point;
>
>
> std::cout << "sphere pull_back: " << std::endl;
> std::cout << q[0] << " " << q[1] << " " << q[2] << std::endl;
> std::cout << "r theta phi" << std::endl;
> std::cout << "..........." << std::endl;
>
> chart_point[0] = q[1];
> chart_point[1] = q[2];
>
> // return (theta,phi) in the chart domain
> return chart_point;
>
> }
> template <int dim,int spacedim, int chartdim>
> Point<spacedim> Ellipsoid<dim,spacedim,chartdim>::push_forward(const
> Point<chartdim> &chart_point) const
> {
> double theta,phi, x,y,z;
>
> phi = chart_point[0];
> theta = chart_point[1];
>
>
> Point<spacedim> p(max_axis,theta,phi);
> // map theta,phi in chart domain onto reference_sphere with radius max_axis
> Point<spacedim> X = reference_sphere.push_forward(p);
>
> // map point on sphere onto ellipsoid
>
> x = a*X[0];
> y = b*X[1];
> z = c*X[2];
>
> Point<spacedim> space_point(x,y,z);
>
> // return point on ellipsoid
> return space_point;
> }
>
>
> Point<3> grid_transform (const Point<3> &X)
> {
> // transform points on sphere onto ellipsoid
> double a,b,c;
> a = 1.0; b = 1.0; c = 1.0;
>
> double x,y,z;
> x = a*X(0);
> y = b*X(1);
> z = c*X(2);
>
> return Point<3>(x,y,z);
> }
>
> void assemble_mesh_and_manifold()
> {
>
> const int dim = 2;
> const int spacedim = 3;
> const int chartdim = 2;
>
> double a,b,c;
> a = 1; b=1; c=1;
>
> Ellipsoid<dim,spacedim,chartdim> ellipsoid(a,b,c);
>
> Triangulation<dim,spacedim> tria;
>
> // generate coarse spherical mesh
> GridGenerator::hyper_sphere (tria, Point<spacedim>(0.0,0.0,0.0), 1.0);
> for (Triangulation<dim,spacedim>::active_cell_iterator
> cell=tria.begin_active(); cell!=tria.end(); ++cell)
> cell->set_all_manifold_ids(0);
>
> print_mesh_info(tria, "spherical_mesh.vtk");
>
> GridTools::transform(&grid_transform, tria);
>
> tria.set_manifold(0,ellipsoid);
>
> tria.refine_global(1);
>
> print_mesh_info(tria, "ellipsoidal_mesh.vtk");
>
> }
>
> int main ()
> {
> assemble_mesh_and_manifold();
> }
>
> Attached is the .vtk of my refined ellipsoidal mesh and the .cc file that
> generates this output (mostly reproduced above.
>
> 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 dealii+unsubscr...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <sphere_to_ellipsoid.cc><ellipsoidal_mesh.vtk>
--
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
