Hello
I'd like to split a triangulation into two and save one of them to file as
a ucd mesh. At present, this is what I have, which throws up an error I'm
not sure how to fix:
std::vector< Point<3> > new_vertices;
std::vector< unsigned int > local_dof_indices( fe.dofs_per_cell );
std::vector< std::vector<unsigned int> > cell_vertices;
std::vector< unsigned int > cell_material_id;
Triangulation<3> new_triangulation;
unsigned int material;
// Loop through all cells
int count = 0;
typename DoFHandler<3>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for ( ; cell != endc; ++cell )
{
// Obtain material id
material = cell->material_id();
if ( (material >= 1) and (material <= 12) ) {
++count;
cell->get_dof_indices( local_dof_indices );
cell_vertices.push_back( local_dof_indices );
cell_material_id.push_back( material );
for (unsigned int v = 0; v <
GeometryInfo<3>::vertices_per_cell; ++v) {
Point<3> vertex_coord( cell->vertex(v) );
new_vertices.push_back( vertex_coord );
}
}
}
std::cout << "Num new cells: " << cell_vertices.size() << std::endl;
std::cout << "Num new cells: " << count << std::endl;
std::vector< CellData<3> > new_cells( cell_vertices.size(),
CellData<3>() );
for ( unsigned int i = 0; i < cell_vertices.size(); ++i ) {
for ( unsigned int j = 0; j <
GeometryInfo<3>::vertices_per_cell; ++j ) {
new_cells[i].vertices[j] = cell_vertices[i][j];
}
new_cells[i].material_id = cell_material_id[i];
}
new_triangulation.create_triangulation( new_vertices, new_cells,
SubCellData() );
This gives the error:
----------------------------------------------------
Exception on processing:
--------------------------------------------------------
An error occurred in line <1928> of file
</home/tedkord/deal.II/source/grid/tria.cc> in function
static void
dealii::internal::Triangulation::Implementation::create_triangulation(const
std::vector<dealii::Point<dim, double> >&, const
std::vector<dealii::CellData<3>, std::allocator<dealii::CellData<3> > >&,
const dealii::SubCellData&, dealii::Triangulation<3, spacedim>&) [with int
spacedim = 3]
The violated condition was:
dealii::GridTools::cell_measure(triangulation.vertices,
cells[cell_no].vertices) >= 0
The name and call sequence of the exception was:
ExcGridHasInvalidCell(cell_no)
Additional Information:
Something went wrong when making cell 2. Read the docs and the source code
for more information.
Many thanks
--
TK
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
