Dear all I am trying to use OpenCASCADE::NormalToMeshProjectionManifold for a STEP file and an initial mesh for an ellipsoid x^2 + 16y^2 + 64z^2 = 1.
I make the step file in GMSH and I use an appropriately scaled unrefined
hyper_sphere as the initial mesh. The error I am getting is
An error occurred in line <241> of file
<./source/opencascade/manifold_lib.cc> in function
dealii::Point<dim>
dealii::OpenCASCADE::{anonymous}::internal_project_to_manifold(const
TopoDS_Shape&, double, const dealii::ArrayView<const dealii::Point<dim> >&,
const dealii::Point<
dim>&) [with int spacedim = 3]
The violated condition was:
closest_point(sh, surrounding_points[i], tolerance)
.distance(surrounding_points[i]) < std::max(tolerance *
surrounding_points[i].norm(), tolerance)
Additional information:
The point [ -0.57735 -0.144338 -0.0721688 ] is not on the manifold.
If I substitute the above point coordinates in the equation of my surface,
I get
x^2 + 16y^2 + 64z^2 = 1.00000183878416
So the point should be on the manifold upto tolerance. I have tried
reducing the tolerance when reading OpenCASCADE::get_shape_tolerance() but
it has not helped.
Please help me debug this.
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/9d6e9a89-f7df-4711-9ea2-5ef39e411fcen%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")
set(TARGET "cad_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_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/opencascade/manifold_lib.h>
#include <deal.II/opencascade/utilities.h>
#include <cmath>
using namespace dealii;
int main(int argc, char **argv) {
// Read the CAD geometry
// TopoDS_Shape ellipsoid = OpenCASCADE::read_IGES("../ellipsoid.iges");
TopoDS_Shape ellipsoid = OpenCASCADE::read_STEP("ellipsoid.step");
const double tolerance = OpenCASCADE::get_shape_tolerance(ellipsoid) * 100;
std::cout << "Shape tolerance = " << tolerance << std::endl;
std::vector<TopoDS_Compound> compounds;
std::vector<TopoDS_CompSolid> compsolids;
std::vector<TopoDS_Solid> solids;
std::vector<TopoDS_Shell> shells;
std::vector<TopoDS_Wire> wires;
OpenCASCADE::extract_compound_shapes(ellipsoid, compounds, compsolids, solids,
shells, wires);
std::cout << "Compounds:" << compounds.size();
std::cout << "\nCompSolids:" << compsolids.size();
std::cout << "\nSolids:" << solids.size();
std::cout << "\nShells:" << shells.size();
std::cout << "\nWires:" << wires.size() << std::endl;
// Define new manifold ids for the triangulation
OpenCASCADE::NormalToMeshProjectionManifold<2, 3> cell_projector(shells[0],
tolerance);
//OpenCASCADE::ArclengthProjectionLineManifold<2, 3> line_projector(wires[0],
// tolerance);
// Make the manifold
double R = 1;
double A = 1;
double B = 16;
double C = 64;
// 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 manifold on the unrefined hypersphere
for (auto &cell : triangulation.active_cell_iterators()) {
cell->set_manifold_id(0);
for (const auto &face : cell->face_iterators()) {
//face->set_manifold_id(1);
face->set_manifold_id(0);
}
}
triangulation.set_manifold(0, cell_projector);
//triangulation.set_manifold(1, line_projector);
// Refine the mesh and see what happens
for (int i = 0; i < 6; ++i) {
// Refine the mesh globally
if (i > 0) {
triangulation.refine_global(i);
}
// Output the final mesh
const std::string filename = "mesh_" + std::to_string(i) + ".vtu";
std::ofstream outfile(filename);
GridOut grid_out;
grid_out.write_vtu(triangulation, outfile);
}
return 0;
}
ellipsoid.step
Description: model/step
