Hi, all.
I'm trying to create a 3D mesh from CAD file (stanford bunny) to solve a
heating equation in this 3D volume.
I tried gmsh, but failed to produce coarse mesh, every try resulted in too
fine meshes (700k cells and more).
Than i tried manually create initial-mesh with one cell (because
`OpenCASCADE::create_triangulation` fails on 3D case). I used FreeCAD to
examine .stp file, selected 8 edges and obtained its vertex coordinates
using python console, and than used a code from Step14
<https://www.dealii.org/current/doxygen/deal.II/step_14.html> Exercise_2_3
to create triangulation.
Finally, i failed to set manifolds correctly.
I have found similar discussions Bruno & Luca
<https://groups.google.com/g/dealii/c/1Y5-iccdhMM/m/R4vkAhF5AQAJ>, yy.wayne
& Wolfgang <https://groups.google.com/g/dealii/c/jgCSqRU550k/m/CVsu7_TlAQAJ>
and
wiki page
<https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-ensure-that-my-refined-grid-is-correct-when-using-cad-geometry>
but
failed to understand the way i should modify sources (see attachment, it is
slightly modified step-54).
Thanks in advance!
Best, Andrey.
--
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 on the web visit
https://groups.google.com/d/msgid/dealii/46a757a1-33be-42a6-898c-d025246564c5n%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2009 - 2021 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
* Authors: Andrea Mola, Luca Heltai, 2014
*/
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/opencascade/manifold_lib.h>
#include <deal.II/opencascade/utilities.h>
// Finally, a few C++ standard header files
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
namespace Step54
{
using namespace dealii;
class TriangulationOnCAD
{
public:
enum ProjectionType
{
NormalProjection = 0,
DirectionalProjection = 1,
NormalToMeshProjection = 2
};
TriangulationOnCAD(
const std::string & initial_mesh_filename,
const std::string & cad_file_name,
const std::string & output_filename,
const ProjectionType surface_projection_kind = NormalProjection);
void run();
private:
void read_domain();
void refine_mesh();
void output_results(const unsigned int cycle);
void create_coarse_grid (Triangulation<3, 3> & coarse_grid);
Triangulation<3, 3> tria;
const std::string initial_mesh_filename;
const std::string cad_file_name;
const std::string output_filename;
const ProjectionType surface_projection_kind;
};
TriangulationOnCAD::TriangulationOnCAD(
const std::string & initial_mesh_filename,
const std::string & cad_file_name,
const std::string & output_filename,
const ProjectionType surface_projection_kind)
: initial_mesh_filename(initial_mesh_filename)
, cad_file_name(cad_file_name)
, output_filename(output_filename)
, surface_projection_kind(surface_projection_kind)
{}
void TriangulationOnCAD::read_domain()
{
TopoDS_Shape bow_surface = OpenCASCADE::read_STEP(cad_file_name, 1e-3);
const double tolerance = OpenCASCADE::get_shape_tolerance(bow_surface) * 10;
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(
bow_surface, compounds, compsolids, solids, shells, wires);
this->create_coarse_grid(this->tria);
output_results(0);
for (const auto &cell : tria.active_cell_iterators())
{
cell->set_all_manifold_ids(1);
for (const auto &face : cell->face_iterators())
{
face->set_all_manifold_ids(2);
for (auto k: face->line_indices())
{
face->line(k)->set_all_manifold_ids(1);
}
}
}
Assert(
wires.size() > 0,
ExcMessage(
"I could not find any wire in the CAD file you gave me. Bailing out."));
OpenCASCADE::ArclengthProjectionLineManifold<3, 3> line_projector(
wires[0], tolerance);
tria.set_manifold(2, line_projector);
{
{
OpenCASCADE::NormalToMeshProjectionManifold<3, 3>
normal_to_mesh_projector(bow_surface, tolerance);
tria.set_manifold(1, normal_to_mesh_projector);
}
}
}
void TriangulationOnCAD::create_coarse_grid (Triangulation<3, 3> & coarse_grid)
{
constexpr auto dim = 3;
const std::vector<Point<dim>> vertices = {
{-17.46951979034 , -16.14238018005 , 0.5628732107468},
{-17.085011860989997 , 3.081818166938 , 0.2961034222146},
{21.51828192439 , -0.3289782258058 , 1.4677078905099998},
{31.57451950505 , -17.329459809190002 , 6.472769511678},
{-28.225385654140002 , 31.51351999628 , 80.75163922927001},
{2.4417120000840002 , 14.50045230315 , 83.30882539870001},
{15.516710343309999 , 16.448021636 , 29.732106105150002},
{28.73074527008 , -10.10172732267 , 45.00440878454}
};
const std::vector<std::array<int, GeometryInfo<dim>::vertices_per_cell>>
cell_vertices = {{{0, 1, 2, 3, 4, 5, 6, 7}}};
auto cells = std::vector<CellData<dim>>(cell_vertices.size(), CellData<dim>{});
for (unsigned int i = 0; i < cell_vertices.size(); ++i)
{
for (unsigned int j = 0; j < cell_vertices[i].size(); ++j)
{
cells[i].vertices[j] = cell_vertices[i][j];
}
cells[i].material_id = 0;
}
coarse_grid.create_triangulation(vertices, cells, SubCellData());
}
void TriangulationOnCAD::refine_mesh()
{
tria.refine_global(1);
}
void TriangulationOnCAD::output_results(const unsigned int cycle)
{
const std::string filename =
(output_filename + "_" + Utilities::int_to_string(cycle) + ".vtk");
std::ofstream logfile(filename);
GridOut grid_out;
grid_out.write_vtk(tria, logfile);
}
void TriangulationOnCAD::run()
{
read_domain();
const unsigned int n_cycles = 5;
for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
{
refine_mesh();
output_results(cycle + 1);
}
}
} // namespace Step54
int main()
{
try
{
using namespace Step54;
const std::string cad_file_name = "input/bunny.stp";
std::string out_mesh_filename = ("3d_mesh_normal_to_mesh_projection");
TriangulationOnCAD tria_on_cad_norm_to_mesh(
"",
cad_file_name,
out_mesh_filename,
TriangulationOnCAD::NormalToMeshProjection);
tria_on_cad_norm_to_mesh.run();
std::cout << "----------------------------------------------------------"
<< std::endl;
std::cout << std::endl;
std::cout << std::endl;
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
