Using CharManifold to implement an Ellipsoid surface gives strange results

[Amit Singh]

Hi

I wrote a code to implement a manifold of an ellipsoid surface embedded in 3d space using 6 chart manifolds corresponding to the 6 faces of a cube. So if the ellipsoid equation is $Ax^2 + By^2 + Cz^2 = R^2$, then I have the following equations:

z_plus(x, y) = sqrt(R^2 - Ax^2 - By^2) / C

z_minus(x, y) = - sqrt(R^2 - Ax^2 - By^2) / C

Similarly, I have analogous equations for x_plus(y, z), x_minus(y, z), y_plus(z, x), and y_plus(z, x).

1. I make an unrefined hyper_sphere<2,3> mesh which is just a cube.

2. Then I transform the sphere into an ellipsoid based on the values of A, B, and C.

3. I reset the SphericalManifold from the triangulation.

4. I detect each of the faces using the face centers and attach one of the six manifolds explained above to each of the faces.

5. Then I refine the mesh globally. The manifold id of the children faces should be inherited from the original 6 faces.

Expected result:

After a few global refinements, I should get something like an ellipsoid.

The result I get:

I have attached my code and the CMakeLists.txt. Please help me understand the error.

Regards,

Amit

The information contained in this electronic communication is intended solely for the individual(s) or entity to which it is addressed. It may contain proprietary, confidential and/or legally privileged information. Any review, retransmission, dissemination, printing, copying or other use of, or taking any action in reliance on the contents of this information by person(s) or entities other than the intended recipient is strictly prohibited and may be unlawful. If you have received this communication in error, please notify us by responding to this email or telephone and immediately and permanently delete all copies of this message and any attachments from your system(s). The contents of this message do not necessarily represent the views or policies of BITS Pilani.

[Amit Singh]

I found my bugs!

There was a mistake in the equation for the push_forward for y_minus chart, and there was also a mistake where I didn't take sqrt of one of the coefficients.

I am able to get Ellipsoidal shapes but the meshes have highly distorted elements at some points as shown below. Any suggestions for a better implementation for this manifold will be helpful. I don't really need the Manifold for mesh refinement but only for higher order mapping.

[Wolfgang Bangerth]

On 1/30/25 08:56, Amit Singh wrote:
>
> I am able to get Ellipsoidal shapes but the meshes have highly distorted
> elements at some points as shown below. Any suggestions for a better
> implementation for this manifold will be helpful.

From the picture, it actually looks like a fairly decent mesh. Can you
point out which points you are unhappy with?

The way I would approach this is to see what happens if you choose
A=B=C, i.e., a sphere. You know how the mesh on a sphere should look
like -- namely, the same way deal.II creates it with the
SphericalManifold. You should then get in essence the same mesh. If you
don't, start looking at why your mesh is different.

Best
W.

[Amit Singh]

Thank you for the reply, @Wolfgang Bangerth. I appreciate your (and the team's) efforts in maintaining deal.ii and even responding to the questions.

I was concerned about some elements of very high aspect ratio occurring at the chart boundaries. They appear like ridges on the surface. I have shown examples for both a sphere and an ellipsoid. Also, the mesh loses symmetry i.e the y-plus and the y-minus patches have different meshes after a few refinements.

Thank you for the tip on debugging.

[Wolfgang Bangerth]

On 1/30/25 22:01, Amit Singh wrote:
>
> I was concerned about some elements of very high aspect ratio occurring at the
> chart boundaries. They appear like ridges on the surface. I have shown
> examples for both a sphere and an ellipsoid. Also, the mesh loses symmetry i.e
> the y-plus and the y-minus patches have different meshes after a few refinements.

I don't know off-hand what is causing this for you, but you probably want to
compare what you are doing with how the GridGenerator::hyper_sphere() function
is doing it:

https://www.dealii.org/current/doxygen/deal.II/namespaceGridGenerator.html#ab2dc2f16fc52cf1e51b9afdc4264ac70
It creates a mesh that does not have this issue, and if you trace where the
difference between your implementation and that of hyper_sphere() comes from,
you will probably also understand how to fix your issue.

Best
W.