obtaining an ellipsoid grid and manifold from GridGenerator::hyper_sphere() and a subclass of ChartManifold

[thomas stephens]

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

[luca.heltai]

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 <tdste...@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 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+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <sphere_to_ellipsoid.cc><ellipsoidal_mesh.vtk>

[thomas stephens]

Luca, Thank you for the help.  

I now have a working Ellipsoid class:

template <int dim,int spacedim>

class Ellipsoid: public SphericalManifold<dim,spacedim>

{

public:

  Ellipsoid(double,double,double);

  Point<spacedim> pull_back(const Point<spacedim> &space_point) const;

  Point<spacedim> push_forward(const Point<spacedim> &chart_point) const;

  Point<spacedim> get_new_point(const Quadrature<spacedim> &quad) const;

  Point<spacedim> grid_transform(const Point<spacedim> &X);

private:

  double  a,b,c;

  double max_axis;

  const Point<spacedim> center;

  Point<dim> ellipsoid_pull_back(const Point<spacedim> &space_point) const;

  Point<spacedim> ellipsoid_push_forward(const Point<dim> &chart_point) const;

};

with member functions: 

template <int dim, int spacedim>

Ellipsoid<dim,spacedim>::Ellipsoid(double a, double b, double c) : SphericalManifold<dim,spacedim>(center), a(a), b(b),c(c), center(0,0,0)

{

  max_axis = std::max(std::max(a,b),c);

}

template <int dim,int spacedim>

Point<spacedim> Ellipsoid<dim,spacedim>::pull_back(const Point<spacedim> &space_point) const

{

  Point<dim> chart_point = ellipsoid_pull_back(space_point);

  Point<spacedim> p;

  p[0] = -1; // dummy radius to match return of SphericalManifold::pull_back()

  p[1] = chart_point[0];

  p[2] = chart_point[1];

  return p;

}

template <int dim,int spacedim>

Point<spacedim> Ellipsoid<dim,spacedim>::push_forward(const Point<spacedim> &chart_point) const

{

  Point<dim> p;  // 

  p[0] = chart_point[1];

  p[1] = chart_point[2];

  Point<spacedim> space_point = ellipsoid_push_forward(p);

  return space_point;

}

template <int dim,int spacedim>

Point<spacedim> Ellipsoid<dim,spacedim>::get_new_point(const Quadrature<spacedim> &quad) const

{

  double u,v,w;

  std::vector< Point<spacedim> > space_points;

  for (unsigned int i=0; i<quad.size(); ++i)

  {

    u = quad.point(i)[0]/a;

    v = quad.point(i)[1]/b;

    w = quad.point(i)[2]/c;

    space_points.push_back(Point<spacedim>(u,v,w));

  }

  Quadrature<spacedim> spherical_quad = Quadrature<spacedim>(space_points, quad.get_weights());

  Point<spacedim> p = SphericalManifold<dim,spacedim>::get_new_point(spherical_quad);

  double x,y,z;

  x = a*p[0];

  y = b*p[1];

  z = c*p[2];

  Point<spacedim> new_point = Point<spacedim>(x,y,z);

  return new_point;

}

template <int dim,int spacedim>

Point<dim> Ellipsoid<dim,spacedim>::ellipsoid_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 = pull_back(p);

  Point<dim> 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>

Point<spacedim> Ellipsoid<dim,spacedim>::ellipsoid_push_forward(const Point<dim> &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 = 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;

}

template<int dim, int spacedim>

Point<spacedim> Ellipsoid<dim,spacedim>::grid_transform(const Point<spacedim> &X)

{

  // transform points X from sphere onto ellipsoid

  double x,y,z;

  x = a*X(0);

  y = b*X(1);

  z = c*X(2);

  return Point<spacedim>(x,y,z);

}

along with the non-member function grid_transform() (which should not be necessary, but I cannot seem to bind the member function Ellipsoid<dim,spacedim>::grid_transform() to an instance of Ellipsoid and a dummy variable, as is done in step-53)

Point<3> grid_transform(const Point<3> &X)

{

  // transform points X from sphere onto ellipsoid

  double x,y,z;

  double a = 1; double b = 3; double c = 5;

  x = a*X(0);

  y = b*X(1);

  z = c*X(2);

  return Point<3>(x,y,z);

}

At any rate, this can all be scripted using the following:

void assemble_mesh_and_manifold()

{

  const int dim = 2;

  const int spacedim = 3;

  double a,b,c;

  a = 1; b=3; c=5;

  Ellipsoid<dim,spacedim> 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);

  //

  //GridTools::transform(std_cxx11::bind(&Ellipsoid<dim,spacedim>::grid_transform,std_cxx11::cref(ellipsoid),std_cxx11::_1), tria); // error when trying to bind to member function in same way as step-53

  tria.set_manifold(0,ellipsoid);

  tria.refine_global(3);

  print_mesh_info(tria, "ellipsoidal_mesh.vtk");

}

I have attached the entire .cc file, and the output from the above code looks great: 

[SebG]

Hi Thomas, hi Luca,

I would like to re-generate the ellipsoidal with deal.ii 9.0.1 but the interfaces of the SphericalManifold class have changed in some release. So I tried to re-implement the following methods

    virtual std::unique_ptr<Manifold<dim,spacedim>> clone() const override;

    virtual Point<spacedim> get_intermediate_point(
            const Point<spacedim> &p1,
            const Point<spacedim> &p2,
            const double weight) const override;

    virtual void get_new_points(
            const ArrayView<const Point<spacedim>>  &surrounding_points,
            const Table<2,double>                   &weights,
            ArrayView<Point<spacedim>>               new_points) const override;

    virtual Point<spacedim> get_new_point(
            const ArrayView<const Point<spacedim>>  &vertices,
            const ArrayView<const double>           &weights) const override;

and two additional methods for the mapping from ellipsoid to sphere:

    /*
     * map a (cartesian) point on a sphere to a point on an ellipsoid
     */
    Point<spacedim>     push_forward(const Point<spacedim>    &space_point) const;

    /*
     * map a (cartesian) point on an ellipsoid to a point on a sphere
     */

    Point<spacedim>     pull_back(const Point<spacedim>  &space_point) const;

Somehow, my implementation does not work well and the mesh is distorted especially for points on the y- and z-axis. I am actually using the ellipsoid as test case and would like generate a mesh consisting of a sphere with (spherical harmonic) boundary topography.

Does anyone know what is going wrong here and the code by Thomas may used in a newer version of deal.ii?

Best wishes,

Sebastian

[luca.heltai]

Did you look at this class?

https://www.dealii.org/developer/doxygen/deal.II/classEllipticalManifold.html

It is in deal.II 9.1.0…

L.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/211cd6b5-3810-4684-89a4-6c43093e125d%40googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.

> <mesh.png><sphere_to_ellipsoid.cc>

[SebG]

Hi Luca,

actually I was not aware of this class. But this class is unfortunately only implemented in 2D. I am not exactly looking for an ellipsoid but for a sphere with a peturbed radius, where r = R * (1 +  e * Y^m_n(\theta, \phi) ). Here R is the reference radius, e > 0 is a scalar and  Y^m_n is a real-valued spherical harmonic.

To implement this, I thought that I could map a point from the perturbed sphere to the reference sphere. Then the updated point is computed on the reference sphere and mapped back to the perturbed one.

As a first step, I tried this approach based on the step-53 tutorial. But I realized that closing the sphere at the poles and in azimuthal direction does not work. Then I thought that it might easier to program the example of an ellipsoid from this discussion.

It would interessting to know why my approach of mapping back to a reference sphere is not working with the updated version of SphericalManifold. But maybe this is worth a new discussion.

Best,

Sebastian

> To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.

[Wolfgang Bangerth]

On 6/11/19 1:06 AM, SebG wrote:
>
> actually I was not aware of this class. But this class is unfortunately only
> implemented in 2D. I am not exactly looking for an ellipsoid but for a sphere
> with a peturbed radius, where r = R * (1 +  e * Y^m_n(\theta, \phi) ). Here R
> is the reference radius, e > 0 is a scalar and  Y^m_n is a real-valued
> spherical harmonic.
>
> To implement this, I thought that I could map a point from the perturbed
> sphere to the reference sphere. Then the updated point is computed on the
> reference sphere and mapped back to the perturbed one.
>
> As a first step, I tried this approach based on the step-53 tutorial. But I
> realized that closing the sphere at the poles and in azimuthal direction does
> not work. Then I thought that it might easier to program the example of an
> ellipsoid from this discussion.
>
> It would interessting to know why my approach of mapping back to a reference
> sphere is not working with the updated version of SphericalManifold. But maybe
> this is worth a new discussion.

I don't know what the issue is, but the way I would write this is a
concatenation of two operations. The first would be the transformation from
reference coordinates to the sphere (done by the push_forward() function of
the SphericalManifold) and the second the scaling in radial direction (done by
your own code). step-53 does this same kind of 2-step procedure, but instead
of using the WGS84 transformation used there, you'd just defer to
SphericalManifold::push_forward.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[SebG]

Hi Wolfgang,

that is exactly the way how I tried to implement this. For a section of a sphere the approach of step-53 works perfectly well for my case too. But I need a full sphere. For the full sphere the problem is that SphericalManifold (like Manifold) does not have the a function called push_forward anymore. So I tried to overload the functions get_intermediate_point, get_new_point and get_new_points. I am not sure which member functions I need to overload if I derive my from class SphericalManifold. I am also not sure if I should derive my class from SphericalManifold or Manifold. It would be great if you could give some advice.





Best wishes,

Sebastian


[luca.heltai]

Dear Sebastian,

have you tried using a simple SphericalManifold for your problem?

The way spherical manifold works is by interpolating *both* the radius and the angles, so that points intermediate between two points with different radii get an average radius between the two.

This may not be the exact solution you were looking for, but if you have a sphere which is uniformly deformed along one axis, then SphericalManifold should do the correct job for you.

L.

> --
> 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+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/5d4bb9c2-d11c-4c4f-aaf9-035352f1c839%40googlegroups.com.

[Wolfgang Bangerth]

On 6/13/19 4:31 AM, SebG wrote:
>
> that is exactly the way how I tried to implement this. For a section of
> a sphere the approach of step-53 works perfectly well for my case too.
> But I need a full sphere. For the full sphere the problem is that
> SphericalManifold (like Manifold) does not have the a function called
> push_forward anymore.

Oh, I see -- SphericalManifold is not derived from the base class that
uses the push_forward/pull_back interface. I didn't remember that...

In that case, Luca's approach is most likely the most useful one.