Dear Wolfgang,
Thanks for looking into it.
I can't guarantee a timeline, but I'll try to take a look at it this
weekend / early next week.
I'll use your code to debug, it is indeed significantly smaller!. I'll try
to see if I can come up with a very small mesh that breaks it down (like
10cells or so) in a similar fashion and try to work through it afterward.
Thanks for the help, i'll keep you posted :)!
Best,
Bruno
On Thursday, 29 November 2018 00:11:11 UTC-5, Wolfgang Bangerth wrote:
>
>
> 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] <javascript:>
> > <mailto:[email protected] <javascript:>>.
> > For more options, visit https://groups.google.com/d/optout.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> <javascript:>
> 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.
