Bruno,
> I have been having some issues when trying to run a "toy" code (a simple
> Steady Navier-Stokes solver) on a case I used to run.
> The boundary group Id 0, which is manually set in my GMSH script
> (vonKarman.geo), is detected perfectly as both a BC and a manifold when I
> export to GMSH 2 and I get the expected results.
> However, if I export the same GMSH file to V4 the same code does not run
> anymore. I guess there is a change in behavior I did not notice in how dealii
> interprets the new format? What am I doing wrong?
> The error I get is (when running in debug mode):
> |
> --------------------------------------------------------
> Anerror occurred inline <10387>of file
> </home/blaisb/codes/dealii/dealii/source/grid/tria.cc>infunction
> voiddealii::Triangulation<dim,spacedim>::set_all_manifold_ids_on_boundary(dealii::types::boundary_id,dealii::types::manifold_id)[withintdim
>
> =2;intspacedim =2;dealii::types::boundary_id
> =unsignedint;dealii::types::manifold_id =unsignedint]
> Theviolated condition was:
> boundary_found
> Additionalinformation:
> Thegiven boundary_id 0isnotdefinedinthisTriangulation!
>
> Stacktrace:
> -----------
> #0 /home/blaisb/codes/dealii/build/lib/libdeal_II.g.so.9.1.0-pre:
> dealii::Triangulation<2, 2>::set_all_manifold_ids_on_boundary(unsigned int,
> unsigned int)
> #1 ./schursteadynavierstokes: SchurNavierStokesSolver<2>::runVonKarman()
> #2 ./schursteadynavierstokes: main
> --------------------------------------------------------
>
> |
I put this piece of code into your source file right after reading in the mesh:
|
for (const auto &cell : triangulation.active_cell_iterators())
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary())
std::cout << "cell=" << cell << ", face=" << f
<< ", boundary_id=" << cell->face(f)->boundary_id()
<< std::endl;
|
The output shows that all boundary faces have boundary ids between 1 and 8.
Of course you already know that from the error message. The question is how
these get there. I have to say that I don't understand the gmsh file format
version 4 -- it's quite complicated with the entities. But do you think you
could try to debug this problem by stepping into the GridIn and seeing where
and why the SubCellData objects are created that assign nonzero boundary_ids
to the boundary faces? The gmsh format 4 reader is new (just a few weeks) and
it would not be a surprise if it still had bugs in it. We would definitely
appreciate any help with this!
A good step would also be to create a smaller mesh that illustrates the issue.
The mesh you have has 400+ cells and 17 entities. It is difficult to read
through the mesh file and understand what each part does.
For all further work, the following, much-reduced program also shows the issue:
|
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_in.h>
int main ()
{
using namespace dealii;
const unsigned int dim = 2;
Triangulation<2> triangulation;
GridIn<dim> grid_in;
grid_in.attach_triangulation (triangulation);
std::ifstream input_file("vonKarman.msh");
grid_in.read_msh(input_file);
for (const auto &cell : triangulation.active_cell_iterators())
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary())
std::cout << "cell=" << cell << ", face=" << f
<< ", boundary_id=" << cell->face(f)->boundary_id()
<< std::endl;
triangulation.set_all_manifold_ids_on_boundary(0,0);
}
|
Best
W.
> However, BC ID 0 as a physical group is defined explicitely in the script. I
> have joined the .h and .cc of my code ( I am sorry it is ugly as hell, it was
> just made for a quick test).
>
> The GMSH file is as follow (also in the tar archive):
> |
>
> y0=0;
> y1=1;
>
> x0 =0.;
> x1 =10;
>
> xc=1;
> yc=0.5;
> r=0.025;
>
> lo =2.0e-1;
> lc =1.0e-1;
>
> Point(0)={x0,y0,0,lo};
> Point(1)={x0,y1,0,lo};
> Point(2)={x1,y0,0,lo};
> Point(3)={x1,y1,0,lo};
>
> Point(4)={xc,yc,0,lc};
> Point(5)={xc-r,yc,0,lc};
> Point(6)={xc,yc+r,0,lc};
> Point(7)={xc+r,yc,0,lc};
> Point(8)={xc,yc-r,0,lc};
>
>
> // Contour Box
> Line(1)={0,1};
> Line(2)={1,3};
> Line(3)={3,2};
> Line(4)={2,0};
>
> // Inner Circle
> Circle(5)={5,4,6};
> Circle(6)={6,4,7};
> Circle(7)={7,4,8};
> Circle(8)={8,4,5};
>
> LineLoop(1)={1,2,3,4};
> LineLoop(2)={5,6,7,8};
>
> // Creates the physical entities
> PhysicalLine(0)={5,6,7,8};
> PhysicalLine(1)={1};
> PhysicalLine(2)={2,4};
> PhysicalLine(3)={3};
>
> PlaneSurface(1)={1,2};
> PhysicalSurface(2)={1};
> RecombineSurface{1};
>
> |
>
>
> The DEALII lines regarding the manifold are :
> |
> grid_in.read_msh(input_file);
> Point<dim,double>circleCenter(1.0,0.5);
> staticconstSphericalManifold<dim>boundary(circleCenter);
> triangulation.set_all_manifold_ids_on_boundary(0,0);
> triangulation.set_manifold (0,boundary);
> |
>
> And the lines linked to the BC are :
> |
> emplate <intdim>
> voidSchurNavierStokesSolver<dim>::setup_dofs ()
> {
> system_matrix.clear();
> pressure_mass_matrix.clear();
>
> dof_handler.distribute_dofs(fe);
>
> std::vector<unsignedint>block_component(dim+1,0);
> block_component[dim]=1;
> DoFRenumbering::component_wise (dof_handler,block_component);
> dofs_per_block.resize (2);
> DoFTools::count_dofs_per_block (dof_handler,dofs_per_block,block_component);
> unsignedintdof_u =dofs_per_block[0];
> unsignedintdof_p =dofs_per_block[1];
>
> FEValuesExtractors::Vectorvelocities(0);
> {
> nonzero_constraints.clear();
> DoFTools::make_hanging_node_constraints(dof_handler,nonzero_constraints);
> VectorTools::interpolate_boundary_values(dof_handler,
> 0,
> ZeroFunction<dim>(dim+1),
> nonzero_constraints,
>
> fe.component_mask(velocities));
>
> if(simulationCase_==TaylorCouette)
> {
> VectorTools::interpolate_boundary_values(dof_handler,
> 1,
> RotatingWall<dim>(),
> nonzero_constraints,
>
> fe.component_mask(velocities));
> }
> if(simulationCase_==BackwardStep)
> {
> VectorTools::interpolate_boundary_values(dof_handler,
> 1,
> PoiseuilleInlet<dim>(),
> nonzero_constraints,
>
> fe.component_mask(velocities));
> }
>
>
>
> if(simulationCase_==vonKarman)
> {
> std::set<types::boundary_id>no_normal_flux_boundaries;
> no_normal_flux_boundaries.insert (2);
> VectorTools::compute_no_normal_flux_constraints (dof_handler,0,
>
> no_normal_flux_boundaries,
>
> nonzero_constraints
> );
> }
>
> if(simulationCase_==vonKarman)
> {
> VectorTools::interpolate_boundary_values(dof_handler,
> 1,
> ConstantXInlet<dim>(),
> nonzero_constraints,
>
> fe.component_mask(velocities));
> }
> }
> nonzero_constraints.close();
>
> {
> zero_constraints.clear();
> DoFTools::make_hanging_node_constraints(dof_handler,zero_constraints);
> VectorTools::interpolate_boundary_values(dof_handler,
> 0,
> ZeroFunction<dim>(dim+1),
> zero_constraints,
>
> fe.component_mask(velocities));
>
> if(simulationCase_==TaylorCouette||simulationCase_==BackwardStep||simulationCase_==vonKarman)
> {
> VectorTools::interpolate_boundary_values(dof_handler,
> 1,
> ZeroFunction<dim>(dim+1),
> zero_constraints,
>
> fe.component_mask(velocities));
> }
>
> if(simulationCase_==vonKarman)
> {
> std::set<types::boundary_id>no_normal_flux_boundaries;
> no_normal_flux_boundaries.insert (2);
> VectorTools::compute_no_normal_flux_constraints (dof_handler,0,
>
> no_normal_flux_boundaries,
> zero_constraints
> );
> }
>
>
> }
> zero_constraints.close();
> std::cout <<" Number of active cells: "
> <<triangulation.n_active_cells()
> <<std::endl
> <<" Number of degrees of freedom: "
> <<dof_handler.n_dofs()
> <<" ("<<dof_u <<'+'<<dof_p <<')'
> <<std::endl;
>
> }
>
> |
>
>
> --
> 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 [email protected]
> <mailto:[email protected]>.
> For more options, visit https://groups.google.com/d/optout.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: [email protected]
www: http://www.math.colostate.edu/~bangerth/
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
