Hi Wolfgang,
I created a simple example but it worked correctly and is attached (I
modified a test code given earlier my Jean-Paul). The difference between
this simple code and my code stated earlier is that between the two checks
of hanging nodes, the interpolate_boundary_values function is not called in
hp_check.cpp. However, in the piece of code sent earlier, there is this
line:
VectorTools::interpolate_boundary_values(dof_handler,
42,
boundary_v,
boundary_values);
This is somehow causing the problem. If I comment this one, the hanging
constraints are correctly. I am trying to put this as well in my simplified
example. In the meantime, may you be you have an idea why this is messing
up. The online documentation is down, so could not check details of the
function.
Best regards
Deepak
On Tue, Sep 27, 2016 at 4:37 PM, Wolfgang Bangerth <[email protected]>
wrote:
> On 09/27/2016 06:44 AM, Deepak Gupta wrote:
>
>> Dear All,
>>
>> Below is a simple piece of code where one finite element has a different
>> p-order compared to the rest. Thus, I expect certain
>> hanging_node_constraints
>> (which is 4) in the output. What I cannot figure out is why I get *zero
>> *constraints when I check the hanging constraints a second time? More can
>> be
>>
>> understood from the code below which forms part of the setup_system()
>> function
>> in my code. The associated output follows the code.
>>
>> boundary_values.clear();
>> dof_handler.distribute_dofs(fe_collection);
>> hanging_node_constraints.clear();
>> DoFTools::make_hanging_node_constraints(dof_handler,
>> hanging_node_constraints);
>> boundary_info(); //created by me for updating rhe boundary
>> indicators
>> //Applying the boundary conditions
>> BoundaryValues<dim> boundary_v;
>> VectorTools::interpolate_boundary_values(dof_handler,
>> 42,
>> boundary_v,
>> boundary_values);
>> //hanging_node_constraints.close();
>> std::cout<< "No.of degrees of freedom: " << dof_handler.n_dofs() <<
>> "\n";
>> std::cout<<"No. of hanging node constraints :
>> "<<hanging_node_constraints.n_constraints()<<std::endl;
>>
>>
>>
>> boundary_values.clear();
>> hanging_node_constraints.clear();
>> DoFTools::make_hanging_node_constraints(dof_handler,
>> hanging_node_constraints);
>> boundary_info(); //created by me for updating rhe boundary indicators
>> //Applying the boundary conditions
>> VectorTools::interpolate_boundary_values(dof_handler,
>> 42,
>> boundary_v,
>> boundary_values);
>> hanging_node_constraints.close();
>> std::cout<< "No.of degrees of freedom: " << dof_handler.n_dofs() <<
>> "\n";
>> std::cout<<"No. of hanging node constraints :
>> "<<hanging_node_constraints.n_constraints()<<std::endl;
>>
>> *OUTPUT
>> No.of degrees of freedom: 6652
>> No. of hanging node constraints : 4
>> No.of degrees of freedom: 6652
>> No. of hanging node constraints : 0
>>
>
> Deepak, this looks wrong but since this is such a well tested part of the
> library, I can only assume that your code is doing something wrong. But
> then, I've been known to suspect the bug somewhere outside the library :-)
> Can you create a small, self-contained testcase that shows the problem?
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> 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/fo
> rum/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.
>
--
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.
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/numerics/data_out.h>
#include <iostream>
#include <fstream>
using namespace dealii;
int main()
{
const int dim = 2;
Triangulation<dim> tria;
const hp::FECollection<dim> fe_collection (FE_Q<dim>(1), FE_Q<dim>(2));
hp::DoFHandler<dim> dof_handler(tria);
// Two cell domain
std::vector<unsigned int> divisions (dim,1);
divisions[0] = 2;
Point<dim> pt_1, pt_2;
for (unsigned int d=0; d<dim; ++d)
pt_2[d] = 1.0;
GridGenerator::subdivided_hyper_rectangle(tria, divisions, pt_1, pt_2);
// Set active FE index (p-refinement level; default index is 0)
dof_handler.begin_active()->set_active_fe_index(1);
dof_handler.distribute_dofs(fe_collection);
// Build constraints
ConstraintMatrix constraints;
constraints.clear();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
constraints.close();
// Print to screen
std::cout
<< "Number of cells: " << tria.n_active_cells() << "\n"
<< "Number of degrees of freedom: " << dof_handler.n_dofs() << "\n"
<< "Number of constraints: " << constraints.n_constraints() << "\n";
//Second round of checking hanging nodes
dof_handler.distribute_dofs(fe_collection);
// Build constraints
constraints.clear();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
constraints.close();
// Print to screen
std::cout
<< "Number of cells: " << tria.n_active_cells() << "\n"
<< "Number of degrees of freedom: " << dof_handler.n_dofs() << "\n"
<< "Number of constraints: " << constraints.n_constraints() << "\n";
// Output p-refinement level
Vector<double> p_refinement (tria.n_active_cells());
unsigned int cell_no = 0;
for (typename hp::DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active();
cell != dof_handler.end(); ++cell, ++cell_no)
{
p_refinement[cell_no] = cell->active_fe_index();
}
DataOut<dim,hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (p_refinement, "p_refinement_level");
data_out.build_patches ();
std::ostringstream filename;
filename << "solution.vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
return 0;
}
