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:
[image: Screenshot from 2025-01-30 19-30-24.png]
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.
--
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.
To view this discussion visit
https://groups.google.com/d/msgid/dealii/3af95a1d-1593-4cf3-ae01-67861baecaf5n%40googlegroups.com.
# #
# CMake script for the step-3 tutorial program:
# #
# Set the name of the project and target:
# set(TARGET "my_code")
set(TARGET "ellipsoid_mesh")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# file(GLOB_RECURSE TARGET_SRC "source/*.cc")
# file(GLOB_RECURSE TARGET_INC "include/*.h")
# set(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
set(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
cmake_minimum_required(VERSION 3.13.4)
find_package(deal.II 9.4.1
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
if(NOT ${deal.II_FOUND})
message(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
endif()
deal_ii_initialize_cached_variables()
project(${TARGET})
deal_ii_invoke_autopilot()
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/manifold.h>
#include <deal.II/grid/tria.h>
#include <cmath>
#include <memory>
using namespace dealii;
/**
* A manifold of the surface Ax^2 + By^2 + Cz^2 = R^2
*/
class EllipsoidalManifold : public ChartManifold<2, 3, 2> {
public:
enum PatchDirection : int {
x_plus,
x_minus,
y_plus,
y_minus,
z_plus,
z_minus
};
EllipsoidalManifold(double R, double A, double B, double C, PatchDirection pd)
: _R(R), _A(A), _B(B), _C(C), _pd(pd) {}
Point<3> push_forward(const Point<2> &chart_point) const override {
double x1 = chart_point[0];
double x2 = chart_point[1];
double x3;
Point<3> output;
switch (_pd) {
case x_plus:
x3 = _push_forward_common(_B, _C, _A, x1, x2, true);
output[0] = x3;
output[1] = x1;
output[2] = x2;
break;
case x_minus:
x3 = _push_forward_common(_B, _C, _A, x1, x2, false);
output[0] = x3;
output[1] = x1;
output[2] = x2;
break;
case y_plus:
x3 = _push_forward_common(_C, _A, _B, x1, x2, true);
output[0] = x2;
output[1] = x3;
output[2] = x1;
break;
case y_minus:
x3 = _push_forward_common(_B, _C, _A, x1, x2, false);
output[0] = x2;
output[1] = x3;
output[2] = x1;
break;
case z_plus:
x3 = _push_forward_common(_A, _B, _C, x1, x2, true);
output[0] = x1;
output[1] = x2;
output[2] = x3;
break;
case z_minus:
x3 = _push_forward_common(_A, _B, _C, x1, x2, false);
output[0] = x1;
output[1] = x2;
output[2] = x3;
break;
}
std::cout << "Push forward: " << chart_point << " -> " << output << std::endl;
return output;
}
Point<2> pull_back(const Point<3> &space_point) const override {
Point<2> output;
switch (_pd) {
case x_plus:
output[0] = space_point[1];
output[1] = space_point[2];
break;
case x_minus:
output[0] = space_point[1];
output[1] = space_point[2];
break;
case y_plus:
output[0] = space_point[2];
output[1] = space_point[0];
break;
case y_minus:
output[0] = space_point[2];
output[1] = space_point[0];
break;
case z_plus:
output[0] = space_point[0];
output[1] = space_point[1];
break;
case z_minus:
output[0] = space_point[0];
output[1] = space_point[1];
break;
}
std::cout << "Pull back: " << space_point << " -> " << output << std::endl;
return output;
}
std::unique_ptr<Manifold<2, 3>> clone() const override {
return std::make_unique<EllipsoidalManifold>(_R, _A, _B, _C, _pd);
}
private:
double _push_forward_common(double coeff1, double coeff2, double coeff3,
double x1, double x2,
bool positiveDirection) const {
double x3 =
std::sqrt(_R * _R - coeff1 * x1 * x1 - coeff2 * x2 * x2) / coeff3;
if (!positiveDirection) {
x3 *= -1;
}
return x3;
}
double _R, _A, _B, _C;
PatchDirection _pd;
};
int main(int argc, char **argv) {
// Make the manifold
double R = 1;
double A = 1;
double B = 16;
double C = 64;
EllipsoidalManifold x_plus_manifold(R, A, B, C, EllipsoidalManifold::x_plus);
EllipsoidalManifold x_minus_manifold(R, A, B, C,
EllipsoidalManifold::x_minus);
EllipsoidalManifold y_plus_manifold(R, A, B, C, EllipsoidalManifold::y_plus);
EllipsoidalManifold y_minus_manifold(R, A, B, C,
EllipsoidalManifold::y_minus);
EllipsoidalManifold z_plus_manifold(R, A, B, C, EllipsoidalManifold::z_plus);
EllipsoidalManifold z_minus_manifold(R, A, B, C,
EllipsoidalManifold::z_minus);
for (int i = 0; i < 3; ++i) {
// Make a spherical mesh
Triangulation<2, 3> triangulation;
GridGenerator::hyper_sphere<3>(triangulation, Point<3>(), R);
// Transform the sphere to an ellipsoid
GridTools::transform(
[A, B, C](const Point<3> &p) -> Point<3> {
return Point<3>(p[0] / std::sqrt(A), p[1] / std::sqrt(B),
p[2] / std::sqrt(C));
},
triangulation);
// Make all the manifolds flat
triangulation.reset_all_manifolds();
// Set new manifolds on the unrefined hypersphere
for (auto &cell : triangulation.active_cell_iterators()) {
Point<3> center = cell->center();
// A hyper-sphere without any global refinement is a rectangular block
// with 6 faces
if (abs(center[1]) < 1e-6 && abs(center[2]) < 1e-6 && center[0] > 0) {
std::cout << "Found x+ face." << std::endl;
cell->set_all_manifold_ids(0);
} else if (abs(center[1]) < 1e-6 && abs(center[2]) < 1e-6 &&
center[0] < 0) {
std::cout << "Found x- face." << std::endl;
cell->set_all_manifold_ids(1);
} else if (abs(center[2]) < 1e-6 && abs(center[0]) < 1e-6 &&
center[1] > 0) {
std::cout << "Found y+ face." << std::endl;
cell->set_all_manifold_ids(2);
} else if (abs(center[2]) < 1e-6 && abs(center[0]) < 1e-6 &&
center[1] < 0) {
std::cout << "Found y- face." << std::endl;
cell->set_all_manifold_ids(3);
} else if (abs(center[0]) < 1e-6 && abs(center[1]) < 1e-6 &&
center[2] > 0) {
std::cout << "Found z+ face." << std::endl;
cell->set_all_manifold_ids(4);
} else if (abs(center[0]) < 1e-6 && abs(center[1]) < 1e-6 &&
center[2] < 0) {
std::cout << "Found z- face." << std::endl;
cell->set_all_manifold_ids(5);
}
}
// Define new manifold ids for the triangulation
triangulation.set_manifold(0, x_plus_manifold);
triangulation.set_manifold(1, x_minus_manifold);
triangulation.set_manifold(2, y_plus_manifold);
triangulation.set_manifold(3, y_minus_manifold);
triangulation.set_manifold(4, z_plus_manifold);
triangulation.set_manifold(5, z_minus_manifold);
// Refine the mesh globally
if (i > 0){
triangulation.refine_global(i);
}
// Output the final mesh
GridOut grid_out;
const std::string filename = "mesh_" + std::to_string(i) + ".vtu";
std::ofstream outfile(filename);
grid_out.write_vtu(triangulation, outfile);
}
return 0;
}
